split SCL

geosentinel/splitSCL.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+import osgeo.ogr as ogr
+import osgeo.osr as osr
+import numpy as np
+import gdal
+import sys
+
+SCL_COLOR = [
+		(0,0,0,255), #NO_DATA
+		(255,0,0,255), #SATURATED_OR_DEFECTIVE
+                (64,64,64,255), #DARK_AREA_PIXELS
+                (102,51,0,255), #CLOUD_SHADOWS
+                (0,255,0,255), #VEGETATION
+		(255,255,0,255), #NOT_VEGETATED
+                (0,0,255,255), #WATER
+		(128,128,128,255), #UNCLASSIFIED
+                (192,192,192,255), #CLOUD_MEDIUM_PROBABILITY
+		(255,255,255,255), #CLOUD_HIGH_PROBABILITY
+		(102,204,255,255), #THIN_CIRRUS
+		(255,153,255,255), #SNOW
+		 ]
+
+#-------------------------------------------------------------------------------
+def splitSCL(filename, outdir):
+	inputImage = gdal.Open(filename)
+	rows = inputImage.RasterYSize
+	cols = inputImage.RasterXSize
+	geotransform =  inputImage.GetGeoTransform()
+	projection = inputImage.GetProjection()
+
+	sr = osr.SpatialReference()
+	sr.ImportFromWkt(projection)
+	driverTiff = gdal.GetDriverByName('GTiff')
+
+	for c in range(len(SCL_COLOR)):
+		outputfile=outdir+filename[:-4]+"_"+str(c)+".tif"
+		output = driverTiff.Create(outputfile, cols, rows, 4, gdal.GDT_Byte)
+		classMask  = inputImage.GetRasterBand(1).ReadAsArray(0,0,cols,rows)==c
+
+		for b in range(1, output.RasterCount+1):
+			output.GetRasterBand(b).SetNoDataValue(0)
+			output.GetRasterBand(b).WriteArray(classMask*SCL_COLOR[c][b-1])
+			output.GetRasterBand(b).FlushCache()
+
+		output.SetGeoTransform(geotransform)
+		output.SetProjection(projection)
+		output = None
+	inputImage = None
+
+
+#-------------------------------------------------------------------------------
+def main(argv):
+	splitSCL(argv[1], argv[2])
+
+if __name__ == "__main__":
+	main(sys.argv)
+

tools/L2ASCL_AreaProcessing.sh

@@ -230,7 +230,26 @@ else
 	mergeL2ASCL_JSON.py $USERDIR"out/mask/SCL/json/" > $USERDIR"scl_data.json"
 fi
 
-#8.-Create Tiles
+#8.- Split SCL
+echo ${RED}"Spliting SCL image..."${NC}
+if [ $BYTILE -ne "0" ]; then
+        echo "BY TILE"
+
+else
+        echo "BY DATE"
+        cd $USERDIR"out/mask/SCL/"
+        if [ ! -d "split" ]; then
+		mkdir split
+	fi
+
+	
+        for f in $(ls  *.tif) ; do
+		echo $f
+		splitSCL.py $f split/
+	done
+fi
+
+#9.-Create Tiles
 echo ${RED}"Creating Tiles..."${NC}
 if [ $BYTILE -ne "0" ]; then
         echo "BY TILE"
@@ -276,7 +295,7 @@ else
 		cd ..
 	done
 
-        cd $USERDIR"out/mask/SCL/"
+        cd $USERDIR"out/mask/SCL/split/"
         if [ ! -d "tiles" ]; then
                 mkdir tiles
         fi