Skip to content

GeoSentinel
GeoSentinel
Commit a308004e authored by Mario Chirinos Colunga's avatar Mario Chirinos Colunga 💬
Browse files
Options

wkt raster

parent c13bbfdf
Hide whitespace changes
Inline Side-by-side
Showing with 24 additions and 23 deletions
geosentinel/rasterWkt.py
View file @ a308004e
...@@ -82,40 +82,41 @@ def rasterWkt(wkt, inputfile, outputfile): ...@@ -82,40 +82,41 @@ def rasterWkt(wkt, inputfile, outputfile):
driverTiff = gdal.GetDriverByName('GTiff') driverTiff = gdal.GetDriverByName('GTiff')
output = driverTiff.Create(outputfile, cols, rows, 3, gdal.GDT_Byte) output = driverTiff.Create(outputfile, cols, rows, 3, gdal.GDT_Byte)
polygonRaster = driverTiff.Create("polygonTmp.tiff", cols, rows, 1, gdal.GDT_Byte) # polygonRaster = driverTiff.Create("polygonTmp.tiff", cols, rows, 1, gdal.GDT_Byte)
#Set up polygon raster # #Set up polygon raster
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.SetGeoTransform(geotransform)
polygonRaster.SetGeoTransform(geotransform) # polygonRaster.SetGeoTransform(geotransform)
polygonRaster.SetProjection(sr.ExportToWkt()) # polygonRaster.SetProjection(sr.ExportToWkt())
# shapefile # # shapefile
driver = ogr.GetDriverByName("ESRI Shapefile") # driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource("myShape.shp") # data_source = driver.CreateDataSource("myShape.shp")
#Burn Polygon # #Burn Polygon
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)
gdal.RasterizeLayer(polygonRaster, [1], layer, burn_values=[255]) # gdal.RasterizeLayer(polygonRaster, [1], layer, burn_values=[255])
polygonRaster.GetRasterBand(1).FlushCache() # polygonRaster.GetRasterBand(1).FlushCache()
#Draw output image #Draw output image
np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows) # np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
np_polygon = getPolygonArray(intput, wkt)
for b in range(1, output.RasterCount+1): for b in range(1, output.RasterCount+1):
np_input = np_polygon & intput.GetRasterBand(b).ReadAsArray(0,0,cols,rows) np_input = np_polygon & intput.GetRasterBand(b).ReadAsArray(0,0,cols,rows)
output.GetRasterBand(b).SetNoDataValue(-99) output.GetRasterBand(b).SetNoDataValue(-99)
......
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