wkt raster

geosentinel/rasterWkt.py

@@ -25,6 +25,47 @@ def createLayer(wkt, layerName="wkt"):
 	feature = None
 	return layer
 #-------------------------------------------------------------------------------
+def getPolygonArray(image, wkt):
+	rows, cols, geotransform = image.RasterYSize, image.RasterXSize, image.GetGeoTransform()
+
+	sr = osr.SpatialReference()
+	sr.ImportFromWkt(image.GetProjection())
+
+	driverTiff = gdal.GetDriverByName('GTiff')
+	polygonRaster = driverTiff.Create("polygonTmp.tiff", 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.SetProjection(sr.ExportToWkt())
+
+	# shapefile
+	driverShp = ogr.GetDriverByName("ESRI Shapefile")
+	data_source = driverShp.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())
+
+	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()
+
+	np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+
+	data_source = None
+	polygonRaster = None
+
+	return np_polygon
+#-------------------------------------------------------------------------------
 def rasterWkt(wkt, inputfile, outputfile):
 	''' 
 	Draw WKT Polygon into a Image