Python osgeo.ogr.OFTInteger() Examples

The following are code examples for showing how to use osgeo.ogr.OFTInteger(). 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: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 6 votes vote down vote up
def _getOGRFieldTypeName(self, ogrOFTType):
#        return ogr.GetFieldTypeName(ogrOFTType)
        if ogrOFTType == ogr.OFTBinary: 
            return 'ogr.OFTBinary'
        elif ogrOFTType == ogr.OFTDate:
            return 'ogr.OFTDate'
        elif ogrOFTType == ogr.OFTDateTime:
            return 'ogr.OFTDateTime'
        elif ogrOFTType == ogr.OFTInteger:
            return 'ogr.OFTInteger'
        elif ogrOFTType == ogr.OFTReal:
            return 'ogr.OFTReal'
        elif ogrOFTType == ogr.OFTString:
            return 'ogr.OFTString'
        else:
            return 'ogr Field Type constant: ' + str(ogrOFTType)

# ............................................... 
Example 2
Project: lidar   Author: giswqs   File: slicing.py    MIT License 6 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 3
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 4
Project: core   Author: lifemapper   File: partnerData.py    GNU General Public License v3.0 5 votes vote down vote up
def _convertType(self, ogrtype):
        if ogrtype == OFTInteger:
            return 'int'
        elif ogrtype == OFTString:
            return 'str'
        elif ogrtype == OFTReal:
            return 'float'
        else:
            raise LMError('Unknown field type {}'.format(ogrtype))
            
    # ............................................... 
Example 5
Project: core   Author: lifemapper   File: layer2.py    GNU General Public License v3.0 5 votes vote down vote up
def verifyField(self, dlocation, ogrFormat, attrname):
        """
        @summary: Read OGR-accessible data and save the features and 
                     featureAttributes on the Vector object
        @param dlocation: Location of the data
        @param ogrFormat: OGR-supported data format code, available at
                             http://www.gdal.org/ogr/ogr_formats.html
        @return: boolean for success/failure 
        @raise LMError: on failure to read data.
        @note: populateStats calls this
        """
        if attrname is None:
            fieldOk = True
        else:
            fieldOk = False
            if dlocation is not None and os.path.exists(dlocation):
                ogr.RegisterAll()
                drv = ogr.GetDriverByName(ogrFormat)
                try:
                    ds = drv.Open(dlocation)
                except Exception, e:
                    raise LMError(['Invalid datasource' % dlocation, str(e)])
                
                lyrDef = ds.GetLayer(0).GetLayerDefn()
                # Make sure given field exists and is the correct type
                for i in range(lyrDef.GetFieldCount()):
                    fld = lyrDef.GetFieldDefn(i)
                    fldname = fld.GetNameRef()
                    if attrname == fldname:
                        fldtype = fld.GetType()
                        if fldtype in (ogr.OFTInteger, ogr.OFTReal, ogr.OFTBinary):
                            fieldOk = True
                        break 
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: 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 8
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 9
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 10
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 11
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 12
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 13
Project: Python-Geospatial-Development-Third-Edition   Author: PacktPublishing   File: utils.py    MIT License 4 votes vote down vote up
def get_ogr_feature_attribute(attr, feature):
    attr_name = attr.name

    if not feature.IsFieldSet(attr_name):
        return (True, None)

    if attr.type == ogr.OFTInteger:
        value = str(feature.GetFieldAsInteger(attr_name))
    elif attr.type == ogr.OFTIntegerList:
        value = repr(feature.GetFieldAsIntegerList(attr_name))
    elif attr.type == ogr.OFTReal:
        value = feature.GetFieldAsDouble(attr_name)
        value = "%*.*f" % (attr.width, attr.precision, value)
    elif attr.type == ogr.OFTRealList:
        values = feature.GetFieldAsDoubleList(attr_name)
        str_values = []
        for value in values:
            str_values.append("%*.*f" % (attr.width,
                                         attr.precision, value))
        value = repr(str_values)
    elif attr.type == ogr.OFTString:
        value = feature.GetFieldAsString(attr_name)
    elif attr.type == ogr.OFTStringList:
        value = repr(feature.GetFieldAsStringList(attr_name))
    elif attr.type == ogr.OFTDate:
        parts = feature.GetFieldAsDateTime(attr_name)
        year,month,day,hour,minute,second,tzone = parts
        value = "%d,%d,%d,%d" % (year,month,day,tzone)
    elif attr.type == ogr.OFTTime:
        parts = feature.GetFieldAsDateTime(attr_name)
        year,month,day,hour,minute,second,tzone = parts
        value = "%d,%d,%d,%d" % (hour,minute,second,tzone)
    elif attr.type == ogr.OFTDateTime:
        parts = feature.GetFieldAsDateTime(attr_name)
        year,month,day,hour,minute,second,tzone = parts
        value = "%d,%d,%d,%d,%d,%d,%d,%d" % (year,month,day,
                                             hour,minute,
                                             second,tzone)
    else:
        return (False, "Unsupported attribute type: " +
                       str(attr.type))

    return (True, value)

############################################################################# 
Example 14
Project: Python-Geospatial-Development-Third-Edition   Author: PacktPublishing   File: utils.py    MIT License 4 votes vote down vote up
def set_ogr_feature_attribute(attr, value, feature):
    attr_name = attr.name

    if value == None:
        feature.UnsetField(attr_name)
        return

    if attr.type == ogr.OFTInteger:
        feature.SetField(attr_name, int(value))
    elif attr.type == ogr.OFTIntegerList:
        integers = eval(value)
        feature.SetFieldIntegerList(attr_name, integers)
    elif attr.type == ogr.OFTReal:
        feature.SetField(attr_name, float(value))
    elif attr.type == ogr.OFTRealList:
        floats = []
        for s in eval(value):
            floats.append(eval(s))
        feature.SetFieldDoubleList(attr_name, floats)
    elif attr.type == ogr.OFTString:
        feature.SetField(attr_name, value)
    elif attr.type == ogr.OFTStringList:
        strings = []
        for s in eval(value):
            strings.append(s.encode(encoding))
        feature.SetFieldStringList(attr_name, strings)
    elif attr.type == ogr.OFTDate:
        parts = value.split(",")
        year  = int(parts[0])
        month = int(parts[1])
        day   = int(parts[2])
        tzone = int(parts[3])
        feature.SetField(attr_name, year, month, day,
                         0, 0, 0, tzone)
    elif attr.type == ogr.OFTTime:
        parts  = value.split(",")
        hour   = int(parts[0])
        minute = int(parts[1])
        second = int(parts[2])
        tzone  = int(parts[3])
        feature.SetField(attr_name, 0, 0, 0,
                         hour, minute, second, tzone)
    elif attr.type == ogr.OFTDateTime:
        parts = value.split(",")
        year   = int(parts[0])
        month  = int(parts[1])
        day    = int(parts[2])
        hour   = int(parts[3])
        minute = int(parts[4])
        second = int(parts[5])
        tzone  = int(parts[6])
        feature.SetField(attr_mame, year, month, day,
                         hour, minute, second, tzone) 
Example 15
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 16
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 17
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 18
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 19
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?