Python osgeo.ogr.FieldDefn() Examples

The following are code examples for showing how to use osgeo.ogr.FieldDefn(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: lidar   Author: giswqs   File: slicing.py    MIT License 7 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)
    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# convert images in a selected folder to shapefiles 
Example 2
Project: swat-tools   Author: DigitalAgricultureDiscovery   File: process.py    GNU General Public License v3.0 6 votes vote down vote up
def update_hru_shapefile(self, inshape, new_data):
        # open hru1.shp
        src = ogr.Open(inshape, 1)

        # get the hru1 layer
        lyr = src.GetLayer()

        # create new field we wish to append (Runoff or Sediment)
        lyr.CreateField(
            ogr.FieldDefn(str(self.fieldswat_output_type).title(), ogr.OFTReal))

        feature = lyr.GetNextFeature()

        for i in range(0, len(new_data)):
            feature.SetField(str(self.fieldswat_output_type), new_data[i])
            lyr.SetFeature(feature)
            feature = lyr.GetNextFeature() 
Example 3
Project: surfclass   Author: Kortforsyningen   File: vectorize.py    MIT License 6 votes vote down vote up
def _add_fields(self):
        self._class_field_index = {}
        vdefn = self._outlyr.GetLayerDefn()
        for class_id, name in self._classmap.items():
            # has the field been created already?
            field_index = vdefn.GetFieldIndex(name)
            if field_index < 0:
                # Create field
                fd = ogr.FieldDefn(name, ogr.OFTInteger64)
                self._outlyr.CreateField(fd)
                field_index = vdefn.GetFieldIndex(name)
            assert field_index >= 0, "Could not create field %s" % name
            self._class_field_index[class_id] = field_index
            logger.debug("Added field: '%s", name)
        # Field for total count
        field_index = vdefn.GetFieldIndex(self._total_field)
        if field_index < 0:
            # Create field
            fd = ogr.FieldDefn(self._total_field, ogr.OFTInteger64)
            self._outlyr.CreateField(fd)
            self._total_field_index = vdefn.GetFieldIndex(self._total_field)
            field_index = vdefn.GetFieldIndex(self._total_field)
            logger.debug("Added field: '%s", self._total_field)
        assert field_index >= 0, "Could not create field %s" % self._total_field 
Example 4
Project: CensusShapefileMaker   Author: GreenInfo-Network   File: GenerateStateMapperData.py    MIT License 6 votes vote down vote up
def main(self):
        print "    Assigning ACS fields to records in censustracts.shp"

        # 1 - open the shapefile and add attributes to it (unless it exists already)
        source = ogr.Open('censusblocks.shp', 1)
        layer  = source.GetLayer()
        defn   = layer.GetLayerDefn()

        if defn.GetFieldIndex('MHHINC') == -1:
            layer.CreateField( ogr.FieldDefn('MHHINC', ogr.OFTInteger) )

        # 2 - then iterate over all records in it, and populate them from the ACS data we cached during init
        # the tract-ID is the first 11 characters of the GEOID, we can key from that
        feature = layer.GetNextFeature()
        while feature:
            tractid    = feature.GetField("GEOID")[:11]
            attributes = self.tract_attributes[tractid]
            feature.SetField("MHHINC", attributes['MHHINC'] )
            layer.SetFeature(feature)
            feature = layer.GetNextFeature()

        # 999 - done
        source.Destroy() 
Example 5
Project: core   Author: lifemapper   File: createshape.py    GNU General Public License v3.0 5 votes vote down vote up
def _addUserFieldDef(self, newDataset):
        spRef = osr.SpatialReference()
        spRef.ImportFromEPSG(DEFAULT_EPSG)
        maxStrlen = LMFormat.getStrlenForDefaultOGR()
     
        newLyr = newDataset.CreateLayer('points', geom_type=ogr.wkbPoint, srs=spRef)
        if newLyr is None:
            raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                    'Layer creation failed')
        
        for idx, vals in self.op.columnMeta.iteritems():
            if vals is not None:
                fldname = str(vals[OccDataParser.FIELD_NAME_KEY])
                fldtype = vals[OccDataParser.FIELD_TYPE_KEY] 
                fldDef = ogr.FieldDefn(fldname, fldtype)
                if fldtype == ogr.OFTString:
                    fldDef.SetWidth(maxStrlen)
                returnVal = newLyr.CreateField(fldDef)
                if returnVal != 0:
                    raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                         'Failed to create field {}'.format(fldname))
                
        # Add wkt field
        fldDef = ogr.FieldDefn(LM_WKT_FIELD, ogr.OFTString)
        fldDef.SetWidth(maxStrlen)
        returnVal = newLyr.CreateField(fldDef)
        if returnVal != 0:
            raise LmException(JobStatus.IO_OCCURRENCE_SET_WRITE_ERROR, 
                                    'Failed to create field {}'.format(fldname))
        
        return newLyr
            

    # ............................................... 
Example 6
Project: lidar   Author: giswqs   File: filling.py    MIT License 5 votes vote down vote up
def polygonize(img,shp_path):
    # mapping between gdal type and ogr field type
    type_mapping = {gdal.GDT_Byte: ogr.OFTInteger,
                    gdal.GDT_UInt16: ogr.OFTInteger,
                    gdal.GDT_Int16: ogr.OFTInteger,
                    gdal.GDT_UInt32: ogr.OFTInteger,
                    gdal.GDT_Int32: ogr.OFTInteger,
                    gdal.GDT_Float32: ogr.OFTReal,
                    gdal.GDT_Float64: ogr.OFTReal,
                    gdal.GDT_CInt16: ogr.OFTInteger,
                    gdal.GDT_CInt32: ogr.OFTInteger,
                    gdal.GDT_CFloat32: ogr.OFTReal,
                    gdal.GDT_CFloat64: ogr.OFTReal}

    ds = gdal.Open(img)
    prj = ds.GetProjection()
    srcband = ds.GetRasterBand(1)

    dst_layername = "Shape"
    drv = ogr.GetDriverByName("ESRI Shapefile")
    dst_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference(wkt=prj)

    dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs)
    # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType])
    raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32])
    dst_layer.CreateField(raster_field)
    gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None)
    del img, ds, srcband, dst_ds, dst_layer


# extract sinks from dem 
Example 7
Project: HackZurich   Author: RefugeeMatchmaking   File: test_shp.py    MIT License 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 8
Project: bitcoin-arbitrage-trading-bot   Author: mammuth   File: jqvmap.py    MIT License 5 votes vote down vote up
def output_ogr(self, output):
    driver = ogr.GetDriverByName( 'ESRI Shapefile' )
    if os.path.exists( output['file_name'] ):
      driver.DeleteDataSource( output['file_name'] )
    source = driver.CreateDataSource( output['file_name'] )
    layer = source.CreateLayer( self.layer_dfn.GetName(),
                                geom_type = self.layer_dfn.GetGeomType(),
                                srs = self.layer.GetSpatialRef() )

    for field in self.fields:
      fd = ogr.FieldDefn( str(field['name']), field['type'] )
      fd.SetWidth( field['width'] )
      if 'precision' in field:
        fd.SetPrecision( field['precision'] )
      layer.CreateField( fd )

    for geometry in self.geometries:
      if geometry.geom is not None:
        feature = ogr.Feature( feature_def = layer.GetLayerDefn() )
        for index, field in enumerate(self.fields):
          if field['name'] in geometry.properties:
            feature.SetField( index, geometry.properties[field['name']].encode('utf-8') )
          else:
            feature.SetField( index, '' )
        feature.SetGeometryDirectly(
          ogr.CreateGeometryFromWkb(
            shapely.wkb.dumps(
              geometry.geom
            )
          )
        )
        layer.CreateFeature( feature )
        feature.Destroy()

    source.Destroy() 
Example 9
Project: qgis-safecast-plugin   Author: OpenGeoLabs   File: __init__.py    GNU General Public License v2.0 5 votes vote down vote up
def _writeToOgrDataSource(self):
        fileExt = self.storageFormat.lower()
        filePath = os.path.splitext(self._fileName)[0] + '.{}'.format(fileExt)
        writer, msg = QgsVectorFileWriter.writeAsVectorFormat(
            self,
            filePath,
            self._provider.encoding(),
            self._provider.crs(),
            driverName=self.storageFormat
        )
        if writer != QgsVectorFileWriter.NoError:
            raise LoadError(
                self.tr("Unable to create SQLite datasource: {}").format(msg)
            )

        # set datasource to new OGR datasource
        self.setDataSource(filePath, self._layerName, 'ogr')
        self._provider = self.dataProvider()

        # create metadata layer
        if self.metadata and 'table' in self.metadata:
            ds = ogr.Open(filePath, True)
            layer_name = self.metadata['table']
            layer = ds.GetLayerByName(layer_name)
            if layer is None:
                layer = ds.CreateLayer(layer_name, None, ogr.wkbNone)
            layer_defn = layer.GetLayerDefn()
            if 'columns' in self.metadata:
                for key in self.metadata['columns']:
                    field = ogr.FieldDefn(key, ogr.OFTString)
                    layer.CreateField(field)

                feat = ogr.Feature(layer_defn)
                for key, value in list(self.metadata['columns'].items()):
                    feat.SetField(key, value)
                layer.CreateFeature(feat)
                feat = None 
Example 10
Project: dzetsaka   Author: nkarasiak   File: function_vector.py    GNU General Public License v3.0 5 votes vote down vote up
def saveToShape(self, array, srs, outShapeFile):
        # Parse a delimited text file of volcano data and create a shapefile
        # use a dictionary reader so we can access by field name
        # set up the shapefile driver
        outDriver = ogr.GetDriverByName('ESRI Shapefile')

        # create the data source
        if os.path.exists(outShapeFile):
            outDriver.DeleteDataSource(outShapeFile)
        # Remove output shapefile if it already exists

        # options = ['SPATIALITE=YES'])
        ds = outDriver.CreateDataSource(outShapeFile)

        # create the spatial reference, WGS84

        lyrout = ds.CreateLayer('randomSubset', srs)
        fields = [
            array[1].GetFieldDefnRef(i).GetName() for i in range(
                array[1].GetFieldCount())]

        for f in fields:
            field_name = ogr.FieldDefn(f, ogr.OFTString)
            field_name.SetWidth(24)
            lyrout.CreateField(field_name)

        for k in array:
            lyrout.CreateFeature(k)

        # Save and close the data source
        ds = None 
Example 11
Project: dzetsaka   Author: nkarasiak   File: function_vector.py    GNU General Public License v3.0 5 votes vote down vote up
def saveToShape(array, srs, outShapeFile):
    # Parse a delimited text file of volcano data and create a shapefile
    # use a dictionary reader so we can access by field name
    # set up the shapefile driver
    import ogr
    outDriver = ogr.GetDriverByName('ESRI Shapefile')

    # create the data source
    if os.path.exists(outShapeFile):
        outDriver.DeleteDataSource(outShapeFile)
    # Remove output shapefile if it already exists

    # options = ['SPATIALITE=YES'])
    ds = outDriver.CreateDataSource(outShapeFile)

    # create the spatial reference, WGS84

    lyrout = ds.CreateLayer('randomSubset', srs, ogr.wkbPoint)
    fields = [array[1].GetFieldDefnRef(i).GetName()
              for i in range(array[1].GetFieldCount())]

    for i, f in enumerate(fields):
        field_name = ogr.FieldDefn(f, array[1].GetFieldDefnRef(i).GetType())
        field_name.SetWidth(array[1].GetFieldDefnRef(i).GetWidth())
        lyrout.CreateField(field_name)

    for k in array:
        lyrout.CreateFeature(k)

    # Save and close the data source
    ds = None 
Example 12
Project: graph-partition   Author: valiantljk   File: test_shp.py    GNU General Public License v2.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c']  # edgenames
        self.paths = (  [(1.0, 1.0), (2.0, 2.0)],
                        [(2.0, 2.0), (3.0, 3.0)],
                        [(0.9, 0.9), (4.0, 2.0)]
                    )
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 13
Project: cadasta-platform   Author: Cadasta   File: shape.py    GNU Affero General Public License v3.0 5 votes vote down vote up
def create_shp_layer(self, layer_type):
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326)

        if layer_type in shp_types.keys():
            layer = self.shp_datasource.CreateLayer(
                layer_type, srs, geom_type=shp_types[layer_type])
            field = ogr.FieldDefn('id', ogr.OFTString)
            layer.CreateField(field)
            return layer 
Example 14
Project: dask-geomodeling   Author: nens   File: test_geometry.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def create_geojson(abspath, polygons=10, bbox=None, ndim=2, projection="EPSG:4326"):
    """Create random triangle polygons inside bbox"""
    driver = ogr.GetDriverByName(str("GeoJSON"))
    driver.DeleteDataSource(abspath)
    datasource = driver.CreateDataSource(abspath)
    layer = datasource.CreateLayer(
        str("results"), get_sr(projection), geom_type=ogr.wkbPolygon
    )
    field_definition = ogr.FieldDefn(str("name"), ogr.OFTString)
    layer.CreateField(field_definition)
    field_definition = ogr.FieldDefn(str("id"), ogr.OFTInteger)
    layer.CreateField(field_definition)
    layer_definition = layer.GetLayerDefn()

    if np.isscalar(polygons):
        polygons = np.random.random((polygons, 3, ndim))
        bbox_min = np.asarray(bbox[:ndim])
        bbox_max = np.asarray(bbox[-ndim:])
        polygons = polygons * (bbox_max - bbox_min) + bbox_min

    for feature_id, coords in enumerate(polygons):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for coord in coords:
            ring.AddPoint_2D(*coord)
        ring.AddPoint_2D(*coords[0])  # close the ring
        polygon = ogr.Geometry(ogr.wkbPolygon)
        polygon.AddGeometry(ring)
        feature = ogr.Feature(layer_definition)
        feature.SetGeometry(polygon)
        feature.SetField(str("name"), str("test"))
        feature.SetField(str("id"), feature_id + 10)
        layer.CreateFeature(feature)
    layer.SyncToDisk()
    datasource.SyncToDisk()

    return polygons 
Example 15
Project: HistoricalMap   Author: nkarasiak   File: function_historical_map.py    GNU General Public License v2.0 5 votes vote down vote up
def polygonize(self,rasterTemp,outShp):
            
            sourceRaster = gdal.Open(rasterTemp)
            band = sourceRaster.GetRasterBand(1)
            driver = ogr.GetDriverByName("ESRI Shapefile")
            # If shapefile already exist, delete it
            if os.path.exists(outShp):
                driver.DeleteDataSource(outShp)
                
            outDatasource = driver.CreateDataSource(outShp)            
            # get proj from raster            
            srs = osr.SpatialReference()
            srs.ImportFromWkt( sourceRaster.GetProjectionRef() )
            # create layer with proj
            outLayer = outDatasource.CreateLayer(outShp,srs)
            # Add class column (1,2...) to shapefile
      
            newField = ogr.FieldDefn('Class', ogr.OFTInteger)
            outLayer.CreateField(newField)
            
            gdal.Polygonize(band, None,outLayer, 0,[],callback=None)  
            
            outDatasource.Destroy()
            sourceRaster=None
            band=None
                    
            
            ioShpFile = ogr.Open(outShp, update = 1)
            
            
            lyr = ioShpFile.GetLayerByIndex(0)
            
            lyr.ResetReading()    

            for i in lyr:
                lyr.SetFeature(i)
            # if area is less than inMinSize or if it isn't forest, remove polygon 
                if i.GetField('Class')!=1:
                    lyr.DeleteFeature(i.GetFID())        
            ioShpFile.Destroy()
                        
            return outShp 
Example 16
Project: uv-align-distribute   Author: c30ra   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 17
Project: honours_project   Author: JFriel   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 18
Project: honours_project   Author: JFriel   File: test_shp.py    GNU General Public License v3.0 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 19
Project: voteswap   Author: sbuss   File: test_shp.py    MIT License 5 votes vote down vote up
def setUp(self):

        def createlayer(driver):
            lyr = shp.CreateLayer("edges", None, ogr.wkbLineString)
            namedef = ogr.FieldDefn("Name", ogr.OFTString)
            namedef.SetWidth(32)
            lyr.CreateField(namedef)
            return lyr

        drv = ogr.GetDriverByName("ESRI Shapefile")

        testdir = os.path.join(tempfile.gettempdir(), 'shpdir')
        shppath = os.path.join(tempfile.gettempdir(), 'tmpshp.shp')

        self.deletetmp(drv, testdir, shppath)
        os.mkdir(testdir)

        shp = drv.CreateDataSource(shppath)
        lyr = createlayer(shp)
        self.names = ['a', 'b', 'c', 'c']  # edgenames
        self.paths = ([(1.0, 1.0), (2.0, 2.0)],
                      [(2.0, 2.0), (3.0, 3.0)],
                      [(0.9, 0.9), (4.0, 0.9), (4.0, 2.0)])

        self.simplified_names = ['a', 'b', 'c']  # edgenames
        self.simplified_paths = ([(1.0, 1.0), (2.0, 2.0)],
                                 [(2.0, 2.0), (3.0, 3.0)],
                                 [(0.9, 0.9), (4.0, 2.0)])

 
        for path, name in zip(self.paths, self.names):
            feat = ogr.Feature(lyr.GetLayerDefn())
            g = ogr.Geometry(ogr.wkbLineString)
            for p in path:
                g.AddPoint_2D(*p)
            feat.SetGeometry(g)
            feat.SetField("Name", name)
            lyr.CreateFeature(feat)
        self.shppath = shppath
        self.testdir = testdir
        self.drv = drv 
Example 20
Project: pyBRAT   Author: Riverscapes   File: bdws.py    GNU General Public License v3.0 5 votes vote down vote up
def setBratFields(self):
        """
        Set values for fields derived from the input BRAT shapefile.

        :return: None
        """
        if self.bratLyr.FindFieldIndex("totdams", 1) >= 0:
            self.bratLyr.DeleteField(self.bratLyr.FindFieldIndex("totdams", 1))
        if self.bratLyr.FindFieldIndex("totcomp", 1) >= 0:
            self.bratLyr.DeleteField(self.bratLyr.FindFieldIndex("totcomp", 1))
        field = ogr.FieldDefn("totdams", ogr.OFTInteger)
        self.bratLyr.CreateField(field)
        field = ogr.FieldDefn("totcomp", ogr.OFTInteger)
        self.bratLyr.CreateField(field) 
Example 21
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 4 votes vote down vote up
def _createPointShapefile(drv, outpath, spRef, lyrname, lyrDef=None,  
                             fldnames=None, idCol=None, xCol=None, yCol=None, 
                             overwrite=True):
        nameChanges = {}
        dlocation = os.path.join(outpath, lyrname + '.shp')
        if os.path.isfile(dlocation):
            if overwrite:
                drv.DeleteDataSource(dlocation)
            else:
                raise LMError('Layer %s exists, creation failed' % dlocation)
        newDs = drv.CreateDataSource(dlocation)
        if newDs is None:
            raise LMError('Dataset creation failed for %s' % dlocation)
        newLyr = newDs.CreateLayer(lyrname, geom_type=ogr.wkbPoint, srs=spRef)
        if newLyr is None:
            raise LMError('Layer creation failed for %s' % dlocation)
        
        # If LayerDefinition is provided, create and add each field to new layer
        if lyrDef is not None:
            for i in range(lyrDef.GetFieldCount()):
                fldDef = lyrDef.GetFieldDefn(i)
#                 fldName = fldDef.GetNameRef()
                returnVal = newLyr.CreateField(fldDef)
                if returnVal != 0:
                    raise LMError('CreateField failed for \'%s\' in %s' 
                                      % (fldDef.GetNameRef(), dlocation))
        # If layer fields are not yet defined, create from fieldnames
        elif (fldnames is not None and idCol is not None and 
                xCol is not None and yCol is not None):
            # Create field definitions
            fldDefList = []
            for fldname in fldnames:
                if fldname in (xCol, yCol):
                    fldDefList.append(ogr.FieldDefn(fldname, ogr.OFTReal))
                elif fldname == idCol:
                    fldDefList.append(ogr.FieldDefn(fldname, ogr.OFTInteger))
                else:
                    fdef = ogr.FieldDefn(fldname, ogr.OFTString)
                    fldDefList.append(fdef)
            # Add field definitions to new layer
            for fldDef in fldDefList:
                try:
                    returnVal = newLyr.CreateField(fldDef)
                    if returnVal != 0:
                        raise LMError('CreateField failed for \'%s\' in %s' 
                                          % (fldname, dlocation))
                    lyrDef = newLyr.GetLayerDefn()
                    lastIdx = lyrDef.GetFieldCount() - 1
                    newFldName = lyrDef.GetFieldDefn(lastIdx).GetNameRef()
                    oldFldName = fldDef.GetNameRef()
                    if newFldName != oldFldName:
                        nameChanges[oldFldName] = newFldName

                except Exception, e:
                    print str(e) 
Example 22
Project: LSDMappingTools   Author: LSDtopotools   File: LSDMap_VectorTools.py    MIT License 4 votes vote down vote up
def geologic_maps_modify_shapefile(shapefile_name, geol_field = "xx"):

	# The shapefile to be rasterized:
	print('Rasterize ' + shapefile_name)
	#get path and filename seperately
	shapefilefilepath = LSDPT.GetPath(shapefile_name)
	#shapefilename = LSDPT.GetFileNameNoPath(shapefile_name)
	shapefileshortname = LSDPT.GetFilePrefix(shapefile_name)

	# get the new shapefile name
	new_shapefile_name = shapefilefilepath+os.sep+shapefileshortname+"_new.shp"

	# copy the shapefile into the new shapefile--we don't wwant to mess up the original data
	print("The New Shapefile name is: "+new_shapefile_name)
	Copy_Shapefile(shapefile_name,new_shapefile_name)

	# New shapefile is opened for writing.
	dataSource = ogr.Open(new_shapefile_name,1)
	daLayer = dataSource.GetLayer(0)

	# add a new field
	new_field = ogr.FieldDefn("GEOL_CODE", ogr.OFTInteger)
	daLayer.CreateField(new_field)

	# lets see what the layers are
	print("Let me tell you what the names of the fields are after I added one!")
	layerDefinition = daLayer.GetLayerDefn()
	for i in range(layerDefinition.GetFieldCount()):
		print(layerDefinition.GetFieldDefn(i).GetName())


	# Make a key for the bedrock
	geol_dict = dict()
	geol_iterator = 0
	geol_field = geol_field #
	for feature in daLayer:
		GEOL = feature.GetField(geol_field)

		if GEOL not in geol_dict:
			geol_iterator = geol_iterator+1
			print("I found a new rock type, GEOL: "+ str(GEOL)+ " and rock type: " + str(geol_iterator))
			geol_dict[GEOL] = geol_iterator

		# now get the geol code
		this_geol_code = geol_dict[GEOL]
		# set the feature
		feature.SetField("GEOL_CODE", this_geol_code)

		# need to update the layer
		daLayer.SetFeature(feature)

	print("The rocks are: ")
	print(geol_dict)

	print("All done")


	return new_shapefile_name, geol_dict 
Example 23
Project: LSDMappingTools   Author: LSDtopotools   File: LSD_GeologyTools.py    MIT License 4 votes vote down vote up
def geologic_maps_modify_shapefile(shapefile_name, geol_field = "xx"):

    # The shapefile to be rasterized:
    print('Rasterize ' + shapefile_name)
    #get path and filename seperately
    shapefilefilepath = LSDPT.GetPath(shapefile_name)
    #shapefilename = LSDPT.GetFileNameNoPath(shapefile_name)
    shapefileshortname = LSDPT.GetFilePrefix(shapefile_name)

    # get the new shapefile name
    new_shapefile_name = shapefilefilepath+os.sep+shapefileshortname+"_new.shp"

    # copy the shapefile into the new shapefile--we don't wwant to mess up the original data
    print("The New Shapefile name is: "+new_shapefile_name)
    Copy_Shapefile(shapefile_name,new_shapefile_name)

    # New shapefile is opened for writing.
    dataSource = ogr.Open(new_shapefile_name,1)
    daLayer = dataSource.GetLayer(0)

    # add a new field
    new_field = ogr.FieldDefn("GEOL_CODE", ogr.OFTInteger)
    daLayer.CreateField(new_field)

    # lets see what the layers are
    print("Let me tell you what the names of the fields are after I added one!")
    layerDefinition = daLayer.GetLayerDefn()
    for i in range(layerDefinition.GetFieldCount()):
        print(layerDefinition.GetFieldDefn(i).GetName())


    # Make a key for the bedrock
    geol_dict = dict()
    geol_iterator = 0
    for feature in daLayer:
        GEOL = feature.GetField(geol_field)

        if GEOL not in geol_dict:
            geol_iterator = geol_iterator+1
            print("I found a new rock type, GEOL: "+ str(GEOL)+ " and rock type: " + str(geol_iterator))
            geol_dict[GEOL] = geol_iterator

        # now get the geol code
        this_geol_code = geol_dict[GEOL]
        # set the feature
        feature.SetField("GEOL_CODE", this_geol_code)

        # need to update the layer
        daLayer.SetFeature(feature)

    print("The rocks are: ")
    print(geol_dict)

    print("All done")


    return new_shapefile_name, geol_dict 
Example 24
Project: swat-tools   Author: DigitalAgricultureDiscovery   File: process.py    GNU General Public License v3.0 4 votes vote down vote up
def create_new_field_shapefile(self, field_shapefile, output_data):

        # set up the shapefile driver
        driver = ogr.GetDriverByName('ESRI Shapefile')

        # create the data source
        data_source = driver.CreateDataSource(self.results_dir + '/Output')

        # create the layer
        layer = data_source.CreateLayer('Field_Response', None, ogr.wkbPolygon)

        # add fields
        layer.CreateField(ogr.FieldDefn('Shape_Area', ogr.OFTReal))
        layer.CreateField(
            ogr.FieldDefn(str(self.fieldswat_output_type), ogr.OFTReal))

        for i in range(0, len(field_shapefile.shapes())):

            # get the field's shape (coordinates)
            field_shape = field_shapefile.shapes()[i]

            # get the field's shape area (float)
            field_record = field_shapefile.records()[i][0]

            # create feature
            feature = ogr.Feature(layer.GetLayerDefn())

            # set the attributes
            feature.SetField('Shape_Area', field_record)
            feature.SetField(str(self.fieldswat_output_type), output_data[i])

            # create polygon ring
            ring = ogr.Geometry(ogr.wkbLinearRing)
            for point in field_shape.points:
                ring.AddPoint(point[0], point[1])

            # create polygon
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            # set the feature geometry using the wkbpolygon
            feature.SetGeometry(poly)

            # create the feature in the layer
            layer.CreateFeature(feature)

            # free up resources
            feature.Destroy()

        # free up resources
        data_source.Destroy() 
Example 25
Project: CensusShapefileMaker   Author: GreenInfo-Network   File: GenerateStateMapperData.py    MIT License 4 votes vote down vote up
def main(self):
        print "    Assigning ACS fields to records in censustracts.shp"

        # 1 - open the shapefile and add attributes to it (unless it exists already)
        source = ogr.Open('censusblocks.shp', 1)
        layer  = source.GetLayer()
        defn   = layer.GetLayerDefn()

        if defn.GetFieldIndex('HISP') == -1:
            layer.CreateField( ogr.FieldDefn('HISP', ogr.OFTInteger) )
        if defn.GetFieldIndex('TOTPOP') == -1:
            layer.CreateField( ogr.FieldDefn('TOTPOP', ogr.OFTInteger) )
        if defn.GetFieldIndex('WHITE') == -1:
            layer.CreateField( ogr.FieldDefn('WHITE', ogr.OFTInteger) )
        if defn.GetFieldIndex('BLACK') == -1:
            layer.CreateField( ogr.FieldDefn('BLACK', ogr.OFTInteger) )
        if defn.GetFieldIndex('AMERIND') == -1:
            layer.CreateField( ogr.FieldDefn('AMERIND', ogr.OFTInteger) )
        if defn.GetFieldIndex('ASIAN') == -1:
            layer.CreateField( ogr.FieldDefn('ASIAN', ogr.OFTInteger) )
        if defn.GetFieldIndex('HAWPI') == -1:
            layer.CreateField( ogr.FieldDefn('HAWPI', ogr.OFTInteger) )
        if defn.GetFieldIndex('YOUTH') == -1:
            layer.CreateField( ogr.FieldDefn('YOUTH', ogr.OFTInteger) )

        # 2 - then iterate over all records in it, and populate them from the ACS data we cached during init
        feature = layer.GetNextFeature()
        while feature:
            geoid    = feature.GetField("GEOID")
            attributes = self.block_attributes[geoid]
            feature.SetField("HISP",    attributes['HISP'] )
            feature.SetField("TOTPOP",  attributes['TOTPOP'] )
            feature.SetField("WHITE",   attributes['WHITE'] )
            feature.SetField("BLACK",   attributes['BLACK'] )
            feature.SetField("AMERIND", attributes['AMERIND'] )
            feature.SetField("ASIAN",   attributes['ASIAN'] )
            feature.SetField("HAWPI",   attributes['HAWPI'] )
            feature.SetField("YOUTH",   attributes['YOUTH'] )
            layer.SetFeature(feature)
            feature = layer.GetNextFeature()

        # 999 - done
        source.Destroy() 
Example 26
Project: TAGGS   Author: jensdebruijn   File: shapefiles.py    MIT License 4 votes vote down vote up
def merge_shapefile_by_id(infile, outfile, layer_name, ID):
    dirname = os.path.dirname(outfile)
    try:
        os.makedirs(dirname)
    except OSError:
        pass
    inshp = ogr.Open(infile)
    inLayer = inshp.GetLayer()

    driver = ogr.GetDriverByName("ESRI Shapefile")

    outshp = driver.CreateDataSource(outfile)
    outLayer = outshp.CreateLayer(layer_name, geom_type=ogr.wkbPolygon)

    idField = ogr.FieldDefn(ID, ogr.OFTInteger)
    outLayer.CreateField(idField)

    gn_ids = [abs(feature.GetField(ID)) for feature in inLayer]
    duplicate_ids = function.find_duplicates(gn_ids)

    duplicate_ids_features = dd(list)

    outLayerDefn = outLayer.GetLayerDefn()

    for i in range(inLayer.GetFeatureCount()):
        inFeature = inLayer.GetFeature(i)
        outFeature = ogr.Feature(outLayerDefn)

        gn_id = abs(inFeature.GetField(ID))
        if gn_id not in duplicate_ids:

            outFeature.SetField(ID, gn_id)

            geom = inFeature.GetGeometryRef()
            outFeature.SetGeometry(geom)

            outLayer.CreateFeature(outFeature)
            outFeature = None

        else:
            duplicate_ids_features[gn_id].append(i)

    for gn_id, features in duplicate_ids_features.items():
        geom = ogr.Geometry(ogr.wkbPolygon)
        for i in features:
            inFeature = inLayer.GetFeature(i)
            geom = geom.Union(inFeature.GetGeometryRef())

        outFeature = ogr.Feature(outLayerDefn)
        outFeature.SetField(ID, gn_id)
        outFeature.SetGeometry(geom)
        outLayer.CreateFeature(outFeature)
        outFeature = None

    inshp = None
    outshp = None 
Example 27
Project: RAPIDpy   Author: erdc   File: voronoi.py    BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def pointsToVoronoiGridShapefile(lat, lon, vor_shp_path, extent=None):
    """
    Converts points to shapefile grid via voronoi
    """
    voronoi_centroids = _get_voronoi_centroid_array(lat, lon, extent)

    # set-up output polygon shp
    log("Creating output polygon shp {0}"
        .format(os.path.basename(vor_shp_path)))
    if os.path.exists(vor_shp_path):
        os.remove(vor_shp_path)
    drv = ogr.GetDriverByName('ESRI Shapefile')
    outShp = drv.CreateDataSource(vor_shp_path)
    osr_geographic_proj = osr.SpatialReference()
    osr_geographic_proj.ImportFromEPSG(4326)
    layer = outShp.CreateLayer('', osr_geographic_proj, ogr.wkbPolygon)
    layer.CreateField(ogr.FieldDefn('GRID_LAT', ogr.OFTReal))
    layer.CreateField(ogr.FieldDefn('GRID_LON', ogr.OFTReal))
    layerDefn = layer.GetLayerDefn()

    # find nodes surrounding polygon centroid
    # sort nodes in counterclockwise order
    # create polygon perimeter through nodes
    log("Building Voronoi polygons...")
    # compute voronoi
    voronoi_manager = Voronoi(voronoi_centroids)
    voronoi_vertices = voronoi_manager.vertices
    voronoi_regions = voronoi_manager.regions
    for point_id, region_index in enumerate(voronoi_manager.point_region):
        vert_index_list = np.array(voronoi_regions[region_index])
        voronoi_centroid = voronoi_centroids[point_id]
        voronoi_poly_points = _get_voronoi_poly_points(vert_index_list,
                                                       voronoi_vertices,
                                                       voronoi_centroid)
        if len(voronoi_poly_points) == 4:
            poly = ogr.Geometry(ogr.wkbPolygon)
            ring = ogr.Geometry(ogr.wkbLinearRing)
            for node in voronoi_poly_points:
                ring.AddPoint(node[0], node[1])

            # grab first node to close ring
            ring.AddPoint(voronoi_poly_points[0][0], voronoi_poly_points[0][1])

            poly.AddGeometry(ring)
            feat = ogr.Feature(layerDefn)
            feat.SetField('GRID_LON', float(voronoi_centroid[0]))
            feat.SetField('GRID_LAT', float(voronoi_centroid[1]))
            feat.SetGeometry(poly)
            layer.CreateFeature(feat) 
Example 28
Project: PyGeoC   Author: lreis2415   File: vector.py    MIT License 4 votes vote down vote up
def raster2shp(rasterfile, vectorshp, layername=None, fieldname=None,
                   band_num=1, mask='default'):
        """Convert raster to ESRI shapefile"""
        FileClass.remove_files(vectorshp)
        FileClass.check_file_exists(rasterfile)
        # this allows GDAL to throw Python Exceptions
        gdal.UseExceptions()
        src_ds = gdal.Open(rasterfile)
        if src_ds is None:
            print('Unable to open %s' % rasterfile)
            sys.exit(1)
        try:
            srcband = src_ds.GetRasterBand(band_num)
        except RuntimeError as e:
            # for example, try GetRasterBand(10)
            print('Band ( %i ) not found, %s' % (band_num, e))
            sys.exit(1)
        if mask == 'default':
            maskband = srcband.GetMaskBand()
        elif mask is None or mask.upper() == 'NONE':
            maskband = None
        else:
            mask_ds = gdal.Open(mask)
            maskband = mask_ds.GetRasterBand(1)
        #  create output datasource
        if layername is None:
            layername = FileClass.get_core_name_without_suffix(rasterfile)
        drv = ogr_GetDriverByName(str('ESRI Shapefile'))
        dst_ds = drv.CreateDataSource(vectorshp)
        srs = None
        if src_ds.GetProjection() != '':
            srs = osr_SpatialReference()
            srs.ImportFromWkt(src_ds.GetProjection())
        dst_layer = dst_ds.CreateLayer(str(layername), srs=srs)
        if fieldname is None:
            fieldname = layername.upper()
        fd = ogr_FieldDefn(str(fieldname), OFTInteger)
        dst_layer.CreateField(fd)
        dst_field = 0
        result = gdal.Polygonize(srcband, maskband, dst_layer, dst_field,
                                 ['8CONNECTED=8'], callback=None)
        return result 
Example 29
Project: rastercube   Author: terrai   File: shputils.py    MIT License 4 votes vote down vote up
def polygon_to_shapefile(polygons, poly_sr, fname, fields_defs=None,
                         poly_attrs=None):
    """
    Write a set of polygons to a shapefile
    Args:
        polygons: a list of (lat, lng) tuples
        fields_defs: The list of fields for those polygons, as a tuple
                     (name, ogr type) for each field. For example :
                     [('field1', ogr.OFTReal), ('field2', ogr.OFTInteger)]
        poly_attrs: A list of dict containing the attributes for each polygon
                        [{'field1' : 1.0, 'field2': 42},
                         {'field1' : 3.0, 'field2': 60}]
    """
    shp_driver = ogr.GetDriverByName("ESRI Shapefile")
    out_ds = shp_driver.CreateDataSource(fname)
    assert out_ds is not None, "Failed to create temporary %s" % fname
    out_layer = out_ds.CreateLayer(fname, poly_sr, geom_type=ogr.wkbPolygon)

    has_attrs = fields_defs is not None
    if has_attrs:
        attrs_name = []
        for field_name, field_type in fields_defs:
            out_layer.CreateField(ogr.FieldDefn(field_name, field_type))
            attrs_name.append(field_name)

    layer_defn = out_layer.GetLayerDefn()
    for i, poly in enumerate(polygons):
        ring = ogr.Geometry(ogr.wkbLinearRing)
        # gdal uses the (x, y) convention => (lng, lat)
        for point in poly:
            ring.AddPoint(point[1], point[0])
        ring.AddPoint(poly[0][1], poly[0][0])  # re-add the start to close
        p = ogr.Geometry(ogr.wkbPolygon)
        p.AddGeometry(ring)

        out_feature = ogr.Feature(layer_defn)
        out_feature.SetGeometry(p)

        if has_attrs:
            attrs = poly_attrs[i]
            for field_name in attrs_name:
                out_feature.SetField(field_name, attrs[field_name])

        out_layer.CreateFeature(out_feature)

    out_feature.Destroy()
    out_ds.Destroy() 
Example 30
Project: pygeotools   Author: dshean   File: geolib.py    MIT License 4 votes vote down vote up
def geom2shp(geom, out_fn, fields=False):
    """Write out a new shapefile for input geometry
    """
    from pygeotools.lib import timelib
    driverName = "ESRI Shapefile"
    drv = ogr.GetDriverByName(driverName)
    if os.path.exists(out_fn):
        drv.DeleteDataSource(out_fn)
    out_ds = drv.CreateDataSource(out_fn)
    out_lyrname = os.path.splitext(os.path.split(out_fn)[1])[0]
    geom_srs = geom.GetSpatialReference()
    geom_type = geom.GetGeometryType()
    out_lyr = out_ds.CreateLayer(out_lyrname, geom_srs, geom_type)
    if fields:
        field_defn = ogr.FieldDefn("name", ogr.OFTString)
        field_defn.SetWidth(128)
        out_lyr.CreateField(field_defn)
        field_defn = ogr.FieldDefn("path", ogr.OFTString)
        field_defn.SetWidth(254)
        out_lyr.CreateField(field_defn)
        #field_defn = ogr.FieldDefn("date", ogr.OFTString)
        #This allows sorting by date
        field_defn = ogr.FieldDefn("date", ogr.OFTInteger)
        field_defn.SetWidth(32)
        out_lyr.CreateField(field_defn)
        field_defn = ogr.FieldDefn("decyear", ogr.OFTReal)
        field_defn.SetPrecision(8)
        field_defn.SetWidth(64)
        out_lyr.CreateField(field_defn)
    out_feat = ogr.Feature(out_lyr.GetLayerDefn())
    out_feat.SetGeometry(geom)
    if fields:
        #Hack to force output extesion to tif, since out_fn is shp
        out_path = os.path.splitext(out_fn)[0] + '.tif'
        out_feat.SetField("name", os.path.split(out_path)[-1])
        out_feat.SetField("path", out_path)
        #Try to extract a date from input raster fn
        out_feat_date = timelib.fn_getdatetime(out_fn)
        if out_feat_date is not None:
            datestamp = int(out_feat_date.strftime('%Y%m%d'))
            #out_feat_date = int(out_feat_date.strftime('%Y%m%d%H%M'))
            out_feat.SetField("date", datestamp)
            decyear = timelib.dt2decyear(out_feat_date)
            out_feat.SetField("decyear", decyear)
    out_lyr.CreateFeature(out_feat)
    out_ds = None
    #return status? 
Example 31
Project: clustre   Author: niels-anders   File: Clustre.py    GNU General Public License v3.0 4 votes vote down vote up
def writeShapefile(lines, FIDs, IDs, label, seeds, sequence, source, previouslyClassified, L, OE, ON):
    shapeData = ogr.Open(lines, 1)
    layer = shapeData.GetLayer()
    FIDs_sorted = FIDs.copy()
    FIDs_sorted.sort()

    FIDs = IDs[FIDs_sorted.tolist()]
    i = 0   
        
    # required fields in shapefile
    newfields = ['label','order','seed','source', 'length', 'OE', 'ON'] #['id','label','order','seed','source']    
    for new_field in newfields:
        if layer.FindFieldIndex(new_field, 0) == -1:     # check if new fields already exist
            new_field = ogr.FieldDefn(new_field, ogr.OFTReal)   
            layer.CreateField(new_field)
    
    feature = layer.GetNextFeature()
    print '# lines classified with label %d: %d' % (label, len(FIDs))
    while feature:
        F = feature.GetFID() 
        if feature.GetField('label') == None:
            feature.SetField('label', 0)
                
        if previouslyClassified[F] == label:
            # remove earlier classification
            feature.SetField('source', 0)
            feature.SetField('label', 0)
            feature.SetField('order', 0)
            feature.SetField('seed', 0)
            feature.SetField('length', 0)        
            feature.SetField('OE', 0)        
            feature.SetField('ON', 0) 
        if i < len(FIDs):
            if F == FIDs[i]:
                feature.SetField('source', source[FIDs_sorted[i]])
                feature.SetField('label', label)
                feature.SetField('order', sequence[FIDs_sorted[i]])
                feature.SetField('seed', seeds[FIDs_sorted[i]])
                feature.SetField('length', L[FIDs_sorted[i]])        
                feature.SetField('OE', OE[FIDs_sorted[i]])        
                feature.SetField('ON', ON[FIDs_sorted[i]])        
                # Cleanup
                i +=1
        layer.SetFeature(feature)
        feature.Destroy()                   # clean cache
        feature = layer.GetNextFeature()