Commit 82103d9b authored by Mario Chirinos Colunga's avatar Mario Chirinos Colunga 💬

wkt raster

parent 50fa6b80
...@@ -32,7 +32,7 @@ def getPolygonArray(image, wkt): ...@@ -32,7 +32,7 @@ def getPolygonArray(image, wkt):
sr.ImportFromWkt(image.GetProjection()) sr.ImportFromWkt(image.GetProjection())
driverTiff = gdal.GetDriverByName('GTiff') driverTiff = gdal.GetDriverByName('GTiff')
polygonRaster = driverTiff.Create("polygonTmp.tif", 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)
...@@ -79,9 +79,9 @@ def rasterWkt(wkt, inputfile, outputfile): ...@@ -79,9 +79,9 @@ def rasterWkt(wkt, inputfile, outputfile):
sr.ImportFromWkt(inputImage.GetProjection()) sr.ImportFromWkt(inputImage.GetProjection())
# Read the input bands as numpy arrays. # Read the input bands as numpy arrays.
print("driver")
driverTiff = gdal.GetDriverByName('GTiff') driverTiff = gdal.GetDriverByName('GTiff')
output = driverTiff.Create(outputfile, cols, rows, inputImage.RasterCount, gdal.GDT_Byte) output = driverTiff.Create(outputfile, cols, rows, 3, gdal.GDT_Byte)
polygonRaster = driverTiff.Create("polygonTmp.tif", cols, rows, 1, gdal.GDT_Byte) polygonRaster = driverTiff.Create("polygonTmp.tif", cols, rows, 1, gdal.GDT_Byte)
...@@ -94,13 +94,12 @@ def rasterWkt(wkt, inputfile, outputfile): ...@@ -94,13 +94,12 @@ def rasterWkt(wkt, inputfile, outputfile):
polygonRaster.SetProjection(sr.ExportToWkt()) polygonRaster.SetProjection(sr.ExportToWkt())
print("Shapefile")
# shapefile # shapefile
driver = ogr.GetDriverByName("ESRI Shapefile") driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource("myShape.shp") data_source = driver.CreateDataSource("myShape.shp")
print("Burn Polygon")
#Burn Polygon #Burn Polygon
srs = osr.SpatialReference() srs = osr.SpatialReference()
srs.ImportFromEPSG(4326) srs.ImportFromEPSG(4326)
...@@ -108,7 +107,6 @@ def rasterWkt(wkt, inputfile, outputfile): ...@@ -108,7 +107,6 @@ def rasterWkt(wkt, inputfile, outputfile):
geometry = ogr.CreateGeometryFromWkt(wkt) geometry = ogr.CreateGeometryFromWkt(wkt)
feature = ogr.Feature(layer.GetLayerDefn()) feature = ogr.Feature(layer.GetLayerDefn())
print("features")
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)
...@@ -116,12 +114,11 @@ def rasterWkt(wkt, inputfile, outputfile): ...@@ -116,12 +114,11 @@ def rasterWkt(wkt, inputfile, outputfile):
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(inputImage, wkt) # np_polygon = getPolygonArray(inputImage, wkt)
print("Setup")
for b in range(1, output.RasterCount+1): for b in range(1, output.RasterCount+1):
np_input = np_polygon & inputImage.GetRasterBand(b).ReadAsArray(0,0,cols,rows) np_input = inputImage.GetRasterBand(b).ReadAsArray(0,0,cols,rows)
output.GetRasterBand(b).SetNoDataValue(-99) output.GetRasterBand(b).SetNoDataValue(-99)
output.GetRasterBand(b).WriteArray(np_input) output.GetRasterBand(b).WriteArray(np_input)
output.GetRasterBand(b).FlushCache() output.GetRasterBand(b).FlushCache()
...@@ -129,13 +126,13 @@ def rasterWkt(wkt, inputfile, outputfile): ...@@ -129,13 +126,13 @@ def rasterWkt(wkt, inputfile, outputfile):
output.SetGeoTransform(geotransform) output.SetGeoTransform(geotransform)
output.SetProjection(sr.ExportToWkt()) output.SetProjection(sr.ExportToWkt())
print("end")
#Close files #Close files
feature = None feature = None
data_source = None data_source = None
polygonRaster = None
output = None output = None
inputImage = None polygonRaster = None
def main(argv): def main(argv):
rasterWkt(argv[1], argv[2], argv[3]) rasterWkt(argv[1], argv[2], argv[3])
......
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