def createGeoJSONFromRaster(geoJsonFileName, array2d, geom, proj,
layerName="BuildingID",
fieldName="BuildingID"):
memdrv = gdal.GetDriverByName('MEM')
src_ds = memdrv.Create('', array2d.shape[1], array2d.shape[0], 1)
src_ds.SetGeoTransform(geom)
src_ds.SetProjection(proj)
band = src_ds.GetRasterBand(1)
band.WriteArray(array2d)
dst_layername = "BuildingID"
drv = ogr.GetDriverByName("geojson")
dst_ds = drv.CreateDataSource(geoJsonFileName)
dst_layer = dst_ds.CreateLayer(layerName, srs=None)
fd = ogr.FieldDefn(fieldName, ogr.OFTInteger)
dst_layer.CreateField(fd)
dst_field = 1
gdal.Polygonize(band, None, dst_layer, dst_field, [], callback=None)
return
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
Python osgeo.gdal | Polygonize() 实例源码
最新推荐文章于 2023-03-17 08:46:22 发布