update

geosentinel/L2ASCLtoJSON.py

@@ -7,18 +7,30 @@ import gdal
 import sys
 import json
 from collections import Counter
+sys.path.append('/home/mario/git/GeoSentinel')
+from geosentinel import rasterWkt
 #-------------------------------------------------------------------------------
 def L2ASCLtoDict(filename, wkt):
 	''' 
 	L2ASCLtoDict
 	'''
-	intput = gdal.Open(filename)
-	rows, cols, geotransform = intput.RasterYSize, intput.RasterXSize, intput.GetGeoTransform()
-	data = intput.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+	inputImage = gdal.Open(filename)
+	rows, cols, geotransform = inputImage.RasterYSize, inputImage.RasterXSize, inputImage.GetGeoTransform()
+	data = inputImage.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+#	count = Counter([])
+#	for r in data:
+#		count = count + Counter(r)
+#	sclDict = {str(k):v for k,v in dict(count).items()}
+#	print ("data")
+#	print (sclDict)
+
+	data = np.bitwise_and(data, rasterWkt.getPolygonArray(inputImage, wkt))
 	count = Counter([])
 	for r in data:
 		count = count + Counter(r)
 	sclDict = {str(k):v for k,v in dict(count).items()}
+#	print ("data & POLY")
+#	print (sclDict)
 	return sclDict
 #-------------------------------------------------------------------------------
 def main(argv):

geosentinel/rasterWkt.py

@@ -6,86 +6,55 @@ import osgeo.osr as osr
 import numpy as np
 import gdal
 import sys
-#-------------------------------------------------------------------------------
-def createLayer(wkt, layerName="wkt"):
-	# create the spatial reference, WGS84
-	srs = osr.SpatialReference()
-	srs.ImportFromEPSG(4326)
-	layer = data_source.CreateLayer(layerName, srs, ogr.wkbPolygon)
-
-	# Create the point from the Well Known Text
-	geometry = ogr.CreateGeometryFromWkt(wkt)
-	# create the feature
-	feature = ogr.Feature(layer.GetLayerDefn())
-	# Set the feature geometry using the point
-	feature.SetGeometry(geometry)
-	# Create the feature in the layer (shapefile)
-	layer.CreateFeature(feature)
-	# Dereference the feature
-	feature = None
-	return layer
+from collections import Counter
 #-------------------------------------------------------------------------------
 def getPolygonArray(image, wkt_geom):
-	rows, cols, geotransform = image.RasterYSize, image.RasterXSize, image.GetGeoTransform()
-	print(rows, cols, geotransform)
-	sr = osr.SpatialReference(wkt=image.GetProjection())
-#	sr.ImportFromWkt(image.GetProjection())
-
-	# shapefile
-#	driverShp = ogr.GetDriverByName("ESRI Shapefile")
-#	data_source = driverShp.CreateDataSource("myShape.shp")
-
-
-#	srs = osr.SpatialReference()
-#	srs.ImportFromEPSG(4326)
-#	layer = data_source.CreateLayer("wkt", sr, ogr.wkbPolygon)
-#	geometry = ogr.CreateGeometryFromWkt(wkt_geom)
-#	feature = ogr.Feature(layer.GetLayerDefn())
-
-#	feature.SetStyleString("PEN(c:#FF0000,w:5px);")
-#	feature.SetGeometry(geometry)
-#	layer.CreateFeature(feature)
-	
-#	data_source = None
-
-	rast_ogr_ds = ogr.GetDriverByName('Memory').CreateDataSource('out')
-        rast_mem_lyr = rast_ogr_ds.CreateLayer("wkt", sr, ogr.wkbPolygon)
-#	rast_mem_lyr = rast_ogr_ds.CreateLayer( 'poly', srs=sr )
-
-
+	rows = image.RasterYSize
+	cols = image.RasterXSize
+	geotransform =  image.GetGeoTransform()
+	projection = image.GetProjection()
+	sr = osr.SpatialReference(wkt=projection)
+
+	sr4326 = osr.SpatialReference()
+	sr4326.ImportFromEPSG(4326)
+
+
+	driverMEM	= gdal.GetDriverByName('MEM')
+	driverMemory	= ogr.GetDriverByName('Memory')
+	driverTiff	= gdal.GetDriverByName('GTiff')
+        driverShape	= ogr.GetDriverByName("ESRI Shapefile")
+
+	#RASTER
+	target_ds =  driverMEM.Create('raster.tif', cols, rows, 1, gdal.GDT_Byte) #MEM
+	target_ds.GetRasterBand(1).SetNoDataValue(0)
+	target_ds.GetRasterBand(1).WriteArray(np.zeros((rows, cols), dtype=np.uint8))
+	target_ds.SetProjection(projection)
+	target_ds.SetGeoTransform(geotransform)
+	target_ds.GetRasterBand(1).FlushCache()
+
+	#VECTOR
+        rast_ogr_ds = driverMemory.CreateDataSource( 'myshape.shp' )
+	rast_mem_lyr = rast_ogr_ds.CreateLayer( 'POLYGON', sr4326,  ogr.wkbPolygon)
 	feat = ogr.Feature( rast_mem_lyr.GetLayerDefn() )
-	feat.SetGeometry(ogr.CreateGeometryFromWkt(wkt_geom))
-#	feat.SetGeometryDirectly( ogr.Geometry(wkt = wkt_geom) )
-
+	feat.SetGeometry( ogr.CreateGeometryFromWkt(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
-	driverTiff = gdal.GetDriverByName('GTiff')
-	polygonRaster = driverTiff.Create("polygonTmp.tif", cols, rows, 1, gdal.GDT_Byte)
-	polygonRaster.GetRasterBand(1).SetNoDataValue(-99)
-	polygonRaster.GetRasterBand(1).WriteArray(np.zeros((rows, cols)))
-	polygonRaster.GetRasterBand(1).FlushCache()
-	polygonRaster.SetProjection(image.GetProjectionRef())
-	polygonRaster.SetGeoTransform(image.GetGeoTransform()) 
 
+	err = gdal.RasterizeLayer(target_ds, [1], rast_mem_lyr, burn_values=[255])
 
-	#Burn Polygon
-	err = gdal.RasterizeLayer(polygonRaster, [1], rast_mem_lyr, burn_values=[255])
 	if err != 0:
 		print(err)
 		gdaltest.post_reason( 'got non-zero result code from RasterizeLayer' )
 		return 'fail'
 
-	polygonRaster.GetRasterBand(1).FlushCache()
-
-	np_polygon = polygonRaster.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
-
-	rast_ogr_ds = None
-#	data_source = None
-	polygonRaster = None
+	target_ds.GetRasterBand(1).FlushCache()
+	np_polygon = target_ds.GetRasterBand(1).ReadAsArray(0,0,cols,rows)
+	count = Counter([])
+#	for r in np_polygon:
+#		count = count + Counter(r)
+#	colorDict = {str(k):v for k,v in dict(count).items()}
+#	print(colorDict)
+	target_ds = None
 
 	return np_polygon
 #-------------------------------------------------------------------------------
@@ -93,74 +62,32 @@ def rasterWkt(wkt, inputfile, outputfile):
 	''' 
 	Draw WKT Polygon into a Image
 	'''
-	print (wkt)
-#	return
-	print (inputfile)
-	print (outputfile)
+
 	inputImage = gdal.Open(inputfile)
-	rows, cols, geotransform = inputImage.RasterYSize, inputImage.RasterXSize, inputImage.GetGeoTransform()
+	rows = inputImage.RasterYSize
+	cols = inputImage.RasterXSize
+	geotransform =  inputImage.GetGeoTransform()
+	projection = inputImage.GetProjection()
+
 	sr = osr.SpatialReference()
-	sr.ImportFromWkt(inputImage.GetProjection())
-	# Read the input bands as numpy arrays.
+	sr.ImportFromWkt(projection)
 
-	print("driver")
 	driverTiff = gdal.GetDriverByName('GTiff')
 	output = driverTiff.Create(outputfile, cols, rows, inputImage.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())
-
-
-#	print("Shapefile")
-#	# shapefile
-#	driver = ogr.GetDriverByName("ESRI Shapefile")
-#	data_source = driver.CreateDataSource("myShape.shp")
-
 
-#	print("Burn Polygon")
-#	#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())
-
-#	print("features")
-#	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()
-
-##	Draw output image
-#	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 = inputImage.GetRasterBand(b).ReadAsArray(0,0,cols,rows)
-		output.GetRasterBand(b).SetNoDataValue(-99)
-		output.GetRasterBand(b).WriteArray(np_polygon)
+		output.GetRasterBand(b).SetNoDataValue(0)
+		output.GetRasterBand(b).WriteArray(np.bitwise_and(np_polygon, np_input))
 		output.GetRasterBand(b).FlushCache()
 
-#	output.SetGeoTransform(geotransform)
-#	output.SetProjection(sr.ExportToWkt())
-	output.SetProjection(inputImage.GetProjectionRef())
-	output.SetGeoTransform(inputImage.GetGeoTransform()) 
+	output.SetGeoTransform(geotransform)
+	output.SetProjection(projection)
 
-	print("end")
 	#Close files
-#	feature = None
-#	data_source = None
-#	polygonRaster = None
+	polygonRaster = None
 	output = None
 	inputImage = None
 

tools/mergeImages.sh

 #!/bin/sh
-DIR=$1	#Output filename
-POLYGON=$2	#Crop Window
+INDIR=$1	#Output filename
+OUTDIR=$2	#Output filename
+POLYGON=$3	#Crop Window
+DIRNAME=$(echo $INDIR | cut -d"/" -f1)
+echo "INDIR:" $INDIR
+echo "OUTDIR:" $OUTDIR
+
+#echo $POLYGON
+cd $INDIR
+
+#if [ ! -d $OUTDIR ]; then
+#	mkdir $OUTDIR
+#fi
 
 #BOX=$(polygonToBox.py "$POLYGON")
 BOX=$POLYGON
 ##------------------------------------------------------------------------------
-cd $DIR
-DIRNAME=$(echo $DIR | cut -d"/" -f1)
+
+
 echo $DIRNAME
+for f in $(ls);  do echo ${f}; done;
+#exit
 ##------------------------------------------------------------------------------
 MERGEDIMAGE=$DIRNAME"_TCI_60m_merged.tif"
 if [  ! -e $MERGEDIMAGE  ]; then
@@ -17,12 +30,13 @@ else
 	echo PASSING $MERGEDIMAGE FOUND
 fi
 
-CROPEDIMAGE=../merge_out/$DIRNAME"_TCI_60m.tif"
+CROPEDIMAGE=$OUTDIR$DIRNAME"_TCI_60m.tif"
 if [  ! -e $CROPEDIMAGE  ]; then
-	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
+	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of GTiff $MERGEDIMAGE $CROPEDIMAGE
 else
 	echo PASSING $CROPEDIMAGE FOUND
 fi
+
 ##------------------------------------------------------------------------------
 MERGEDIMAGE=$DIRNAME"_SCL_60m_merged.tif"
 if [  ! -e $MERGEDIMAGE  ]; then
@@ -32,9 +46,9 @@ else
 	echo PASSING $MERGEDIMAGE FOUND
 fi
 
-CROPEDIMAGE=../merge_out/$DIRNAME"_SCL_60m.tif"
+CROPEDIMAGE=$OUTDIR$DIRNAME"_SCL_60m.tif"
 if [  ! -e $CROPEDIMAGE  ]; then
-	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
+	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of GTiff $MERGEDIMAGE $CROPEDIMAGE
 else
 	echo PASSING $CROPEDIMAGE FOUND
 fi

tools/mergeImages.sh.save

+#!/bin/sh
+DIR=$1	#Output filename
+POLYGON=$2	#Crop Window
+
+echo $DIR
+echo $POLYGON
+
+
+#BOX=$(polygonToBox.py "$POLYGON")
+BOX=$POLYGON
+##------------------------------------------------------------------------------
+cd $DIR
+DIRNAME=$(echo $DIR | cut -d"/" -f1)
+echo $DIRNAME
+##------------------------------------------------------------------------------
+MERGEDIMAGE=$DIRNAME"_TCI_60m_merged.tif"
+if [  ! -e $MERGEDIMAGE  ]; then
+	echo "gdal_merge.py -o "$MERGEDIMAGE $(ls *TCI_10m.jp2)
+	gdal_merge.py -o $MERGEDIMAGE $(ls *TCI_10m.jp2)
+else
+	echo PASSING $MERGEDIMAGE FOUND
+fi
+
+CROPEDIMAGE=../merge_out/$DIRNAME"_TCI_60m.jp2"
+if [  ! -e $CROPEDIMAGE  ]; then
+	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
+else
+	echo PASSING $CROPEDIMAGE FOUND
+fi
+##------------------------------------------------------------------------------
+MERGEDIMAGE=$DIRNAME"_SCL_60m_merged.tif"
+if [  ! -e $MERGEDIMAGE  ]; then
+	echo "gdal_merge.py -o "$MERGEDIMAGE $(ls *SCL_60m.jp2)
+	gdal_merge.py -o $MERGEDIMAGE $(ls *SCL_60m.jp2)
+else
+	echo PASSING $MERGEDIMAGE FOUND
+fi
+
+CROPEDIMAGE=../merge_out/$DIRNAME"_SCL_60m.jp2"
+if [  ! -e $CROPEDIMAGE  ]; then
+	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
+else
+	echo PASSING $CROPEDIMAGE FOUND
+fi
+
+cd ..
+
+
+

tools/mergeImagesByDirectory.sh

 #!/bin/sh
 JP2DIR=$1	#JP2 Directory
-POLYGON=$2	#Crop Window
-JOBS=${3:-1}
-
-
+OUTDIR=$2
+POLYGON=$3	#Crop Window
+JOBS=${4:-1}
+DIRNAME=$(echo $JP2DIR | cut -d"/" -f1)
+#echo "POLYGON"
+#echo $POLYGON
 BOX=$(polygonToBox.py "$POLYGON")
 if [ "$#" -le 1 ]; then
 	echo "Usage " $0 " <Parent Directory with subdirectories> <POLYGON> [jobs]" 
 	exit
 fi
-
+echo "cd " $JP2DIR
 cd $JP2DIR
-mkdir merge_out
-ls -d */ | parallel -q --jobs $JOBS mergeImages.sh {} "$BOX"
+for f in $(ls);  do echo ${f}; done;
+
+ls -d */ | parallel -q --jobs $JOBS mergeImages.sh {} $OUTDIR "$BOX"
+#ls -d */ | parallel -q --jobs $JOBS test.sh {} ../../{} "$BOX"
+
 
-find . -name *.xml -type f -delete