Commit 5cb902ca authored by Mario Chirinos Colunga's avatar Mario Chirinos Colunga 💬

wkt raster

parent 3569cdd5
...@@ -25,27 +25,40 @@ def createLayer(wkt, layerName="wkt"): ...@@ -25,27 +25,40 @@ def createLayer(wkt, layerName="wkt"):
feature = None feature = None
return layer return layer
#------------------------------------------------------------------------------- #-------------------------------------------------------------------------------
def getPolygonArray(image, wkt): def getPolygonArray(image, wkt_geom):
rows, cols, geotransform = image.RasterYSize, image.RasterXSize, image.GetGeoTransform() rows, cols, geotransform = image.RasterYSize, image.RasterXSize, image.GetGeoTransform()
print(rows, cols, geotransform) print(rows, cols, geotransform)
sr = osr.SpatialReference() sr = osr.SpatialReference()
sr.ImportFromWkt(image.GetProjection()) sr.ImportFromWkt(image.GetProjection())
# shapefile # # shapefile
driverShp = ogr.GetDriverByName("ESRI Shapefile") # driverShp = ogr.GetDriverByName("ESRI Shapefile")
data_source = driverShp.CreateDataSource("myShape.shp") # data_source = driverShp.CreateDataSource("myShape.shp")
srs = osr.SpatialReference() # srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) # srs.ImportFromEPSG(4326)
layer = data_source.CreateLayer("wkt", srs, ogr.wkbPolygon) # layer = data_source.CreateLayer("wkt", srs, ogr.wkbPolygon)
geometry = ogr.CreateGeometryFromWkt(wkt) # geometry = ogr.CreateGeometryFromWkt(wkt)
feature = ogr.Feature(layer.GetLayerDefn()) # feature = ogr.Feature(layer.GetLayerDefn())
# feature.SetStyleString("PEN(c:#FF0000,w:5px);") ## feature.SetStyleString("PEN(c:#FF0000,w:5px);")
feature.SetGeometry(geometry) # feature.SetGeometry(geometry)
layer.CreateFeature(feature) # layer.CreateFeature(feature)
#
# data_source = None
rast_ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource( 'wrk' )
rast_mem_lyr = rast_ogr_ds.CreateLayer( 'poly', srs=sr )
feat = ogr.Feature( rast_mem_lyr.GetLayerDefn() )
feat.SetGeometryDirectly( ogr.Geometry(wkt = wkt_geom) )
rast_mem_lyr.CreateFeature( feat )
# Run the algorithm.
# err = gdal.RasterizeLayer( target_ds, [3,2,1], rast_mem_lyr, burn_values = [200,220,240] )
#Set up polygon raster #Set up polygon raster
driverTiff = gdal.GetDriverByName('GTiff') driverTiff = gdal.GetDriverByName('GTiff')
...@@ -53,18 +66,18 @@ def getPolygonArray(image, wkt): ...@@ -53,18 +66,18 @@ def getPolygonArray(image, wkt):
polygonRaster.GetRasterBand(1).SetNoDataValue(-99) polygonRaster.GetRasterBand(1).SetNoDataValue(-99)
polygonRaster.GetRasterBand(1).WriteArray(np.zeros((rows, cols))) polygonRaster.GetRasterBand(1).WriteArray(np.zeros((rows, cols)))
polygonRaster.GetRasterBand(1).FlushCache() polygonRaster.GetRasterBand(1).FlushCache()
# polygonRaster.SetGeoTransform(geotransform)
# polygonRaster.SetProjection(sr.ExportToWkt())
polygonRaster.SetProjection(image.GetProjectionRef()) polygonRaster.SetProjection(image.GetProjectionRef())
polygonRaster.SetGeoTransform(image.GetGeoTransform()) polygonRaster.SetGeoTransform(image.GetGeoTransform())
#Burn Polygon #Burn Polygon
gdal.RasterizeLayer(polygonRaster, [1], layer, burn_values=[255]) gdal.RasterizeLayer(polygonRaster, [1], rast_mem_lyr, burn_values=[255])
polygonRaster.GetRasterBand(1).FlushCache() polygonRaster.GetRasterBand(1).FlushCache()
np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows) np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
feature = None rast_ogr_ds = None
data_source = None # data_source = None
polygonRaster = None polygonRaster = None
return np_polygon return np_polygon
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment