wkt raster

geosentinel/rasterWkt.py

@@ -73,52 +73,53 @@ def rasterWkt(wkt, inputfile, outputfile):
 	print (wkt)
 	print (inputfile)
 	print (outputfile)
-	intput = gdal.Open(inputfile)
-	rows, cols, geotransform = intput.RasterYSize, intput.RasterXSize, intput.GetGeoTransform()
+	inputImage = gdal.Open(inputfile)
+	rows, cols, geotransform = inputImage.RasterYSize, inputImage.RasterXSize, inputImage.GetGeoTransform()
 	sr = osr.SpatialReference()
-	sr.ImportFromWkt(intput.GetProjection())
+	sr.ImportFromWkt(inputImage.GetProjection())
 	# Read the input bands as numpy arrays.
 
-
+	print("driver")
 	driverTiff = gdal.GetDriverByName('GTiff')
-	output = driverTiff.Create(outputfile, cols, rows, 3, gdal.GDT_Byte)
-	polygonRaster = driverTiff.Create("polygonTmp.tif", cols, rows, 1, gdal.GDT_Byte)
+	output = driverTiff.Create(outputfile, cols, rows, output.RasterCount(), gdal.GDT_Byte)
+#	polygonRaster = driverTiff.Create("polygonTmp.tif", cols, rows, 1, gdal.GDT_Byte)
 
 
-	#Set up polygon raster
-	polygonRaster.GetRasterBand(1).SetNoDataValue(-99)
-	polygonRaster.GetRasterBand(1).WriteArray(np.zeros((rows, cols)))
-	polygonRaster.GetRasterBand(1).FlushCache()
-	polygonRaster.SetGeoTransform(geotransform)
-	polygonRaster.SetGeoTransform(geotransform)
-	polygonRaster.SetProjection(sr.ExportToWkt())
+#	#Set up polygon raster
+#	polygonRaster.GetRasterBand(1).SetNoDataValue(-99)
+#	polygonRaster.GetRasterBand(1).WriteArray(np.zeros((rows, cols)))
+#	polygonRaster.GetRasterBand(1).FlushCache()
+#	polygonRaster.SetGeoTransform(geotransform)
+#	polygonRaster.SetGeoTransform(geotransform)
+#	polygonRaster.SetProjection(sr.ExportToWkt())
 
 
 
-	# shapefile
-	driver = ogr.GetDriverByName("ESRI Shapefile")
-	data_source = driver.CreateDataSource("myShape.shp")
+#	# shapefile
+#	driver = ogr.GetDriverByName("ESRI Shapefile")
+#	data_source = driver.CreateDataSource("myShape.shp")
 
 
-	#Burn Polygon
-	srs = osr.SpatialReference()
-	srs.ImportFromEPSG(4326)
-	layer = data_source.CreateLayer("wkt", srs, ogr.wkbPolygon)
-	geometry = ogr.CreateGeometryFromWkt(wkt)
-	feature = ogr.Feature(layer.GetLayerDefn())
+#	#Burn Polygon
+#	srs = osr.SpatialReference()
+#	srs.ImportFromEPSG(4326)
+#	layer = data_source.CreateLayer("wkt", srs, ogr.wkbPolygon)
+#	geometry = ogr.CreateGeometryFromWkt(wkt)
+#	feature = ogr.Feature(layer.GetLayerDefn())
 
-	feature.SetStyleString("PEN(c:#FF0000,w:5px);")
-	feature.SetGeometry(geometry)
-	layer.CreateFeature(feature)
+#	feature.SetStyleString("PEN(c:#FF0000,w:5px);")
+#	feature.SetGeometry(geometry)
+#	layer.CreateFeature(feature)
 
-	gdal.RasterizeLayer(polygonRaster, [1], layer, burn_values=[255])
-	polygonRaster.GetRasterBand(1).FlushCache()
+#	gdal.RasterizeLayer(polygonRaster, [1], layer, burn_values=[255])
+#	polygonRaster.GetRasterBand(1).FlushCache()
 
 	#Draw output image
-	np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
-#	np_polygon = getPolygonArray(intput, wkt)
+#	np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+#	np_polygon = getPolygonArray(inputImage, wkt)
+	print("Setup")
 	for b in range(1, output.RasterCount+1):
-		np_input = intput.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).WriteArray(np_input)
 		output.GetRasterBand(b).FlushCache()
@@ -126,13 +127,13 @@ def rasterWkt(wkt, inputfile, outputfile):
 	output.SetGeoTransform(geotransform)
 	output.SetProjection(sr.ExportToWkt())
 
-
+	print("end")
 	#Close files
-	feature = None
-	data_source = None
+#	feature = None
+#	data_source = None
+#	polygonRaster = None
 	output = None
-	polygonRaster = None
-
+	inputImage = None
 
 def main(argv):
 	rasterWkt(argv[1], argv[2], argv[3])