merge

geosentinel/L2ASCL_Directory_FillOclusions.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
+#-------------------------------------------------------------------------------
+def L2ASCL_Directory_FillOclusions:
+
+dirList = os.listdir(inDir)
+
+
+#-------------------------------------------------------------------------------
+def main(argv):
+	if len(sys.argv) != 2:
+		print("Usage: " + argv[0] + " <JSON File>")
+	else:
+		jsonFile=open(argv[1]).read()
+		cfg = json.loads(jsonFile)
+		extractSentielData(filename, destination)
+#-------------------------------------------------------------------------------
+if __name__ == "__main__":
+	main(sys.argv)

geosentinel/L2ASCL_FillOclusions.py

@@ -8,11 +8,11 @@ import sys
 
 #-------------------------------------------------------------------------------
 
-def replaceOclussions(filename1, filename2):
+def replaceOclussions(filename1, filename2, outFilename):
 	'''
 		Fill File2 Clouds with File1 Data
 	'''
-	clouds = [3,8,9,10]
+	clouds = [0, 3,8,9,10]
 	input1 = gdal.Open(filename1)
 	rows1, cols1, geotransform1 = input1.RasterYSize, input1.RasterXSize, input1.GetGeoTransform()
 	data1 = input1.GetRasterBand(1).ReadAsArray(0,0,cols1,rows1)
@@ -21,22 +21,35 @@ def replaceOclussions(filename1, filename2):
 	rows2, cols2, geotransform2 = input2.RasterYSize, input2.RasterXSize, input2.GetGeoTransform()
 	data2 = input2.GetRasterBand(1).ReadAsArray(0,0,cols2,rows2)
 
-	for r in range(len(input1)):
-		for c in range(len(input1[r])):
-			if input2[r][c] in clouds:
-				input2[r][c] = input1[r][c]
 
-	input2.GetRasterBand(1).WriteArray(input2)
-	input2.GetRasterBand(1).FlushCache()
+	driverTiff = gdal.GetDriverByName('GTiff')
+	output = driverTiff.Create(outFilename, cols2, rows2, 1, gdal.GDT_Byte)
+
+	sr = osr.SpatialReference()
+	sr.ImportFromWkt(input2.GetProjection())
+	output.GetRasterBand(1).SetNoDataValue(-99)
+	output.GetRasterBand(1).WriteArray(np.zeros((rows2, cols2)))
+	output.GetRasterBand(1).FlushCache()
+	output.SetGeoTransform(geotransform2)
+	output.SetProjection(sr.ExportToWkt())
+
+	for r in range(rows1):
+		for c in range(cols1):
+			if data2.item((r,c)) in clouds:
+				data2.itemset( (r,c), data1.item((r,c)) )
+
+	output.GetRasterBand(1).WriteArray(data2)
+	output.GetRasterBand(1).FlushCache()
 
 	input1 = None
 	input2 = None
+	output = None
+#-------------------------------------------------------------------------------
 def main(argv):
-	if len(sys.argv) != 2:
-		print("Usage: " + argv[0] + " <JSON File>")
+	if len(sys.argv) != 4:
+		print("Usage: " + argv[0] + " <file1> <file2> <output>")
 	else:
-		jsonFile=open(argv[1]).read()
-		cfg = json.loads(jsonFile)
-		extractSentielData(filename, destination)
+		replaceOclussions(argv[1], argv[2], argv[3])
+#-------------------------------------------------------------------------------
 if __name__ == "__main__":
 	main(sys.argv)

geosentinel/out.tif.aux.xml

+<PAMDataset>
+  <PAMRasterBand band="1">
+    <Histograms>
+      <HistItem>
+        <HistMin>-0.4615384615384616</HistMin>
+        <HistMax>12.46153846153846</HistMax>
+        <BucketCount>13</BucketCount>
+        <IncludeOutOfRange>0</IncludeOutOfRange>
+        <Approximate>1</Approximate>
+        <HistCounts>118281|375|1762|6524|28005|9873|8645|11659|1529|6235|569|177|10</HistCounts>
+      </HistItem>
+    </Histograms>
+    <Metadata>
+      <MDI key="STATISTICS_MAXIMUM">12</MDI>
+      <MDI key="STATISTICS_MEAN">2.0369492470719</MDI>
+      <MDI key="STATISTICS_MINIMUM">0</MDI>
+      <MDI key="STATISTICS_STDDEV">2.8042144980274</MDI>
+    </Metadata>
+  </PAMRasterBand>
+</PAMDataset>

tools/mergeImages.sh

@@ -17,7 +17,7 @@ else
 	echo PASSING $MERGEDIMAGE FOUND
 fi
 
-CROPEDIMAGE=../merged/$DIRNAME"_TCI_10m.tif"
+CROPEDIMAGE=../merge_out/$DIRNAME"_TCI_10m.tif"
 if [  ! -e $CROPEDIMAGE  ]; then
 	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
 else
@@ -32,7 +32,7 @@ else
 	echo PASSING $MERGEDIMAGE FOUND
 fi
 
-CROPEDIMAGE=../merged/$DIRNAME"_SCL_20m.tif"
+CROPEDIMAGE=../merge_out/$DIRNAME"_SCL_20m.tif"
 if [  ! -e $CROPEDIMAGE  ]; then
 	gdal_translate -projwin $BOX -projwin_srs WGS84 -ot Byte -of JPEG $MERGEDIMAGE $CROPEDIMAGE
 else

tools/mergeImagesByDirectory.sh

@@ -11,7 +11,7 @@ if [ "$#" -le 1 ]; then
 fi
 
 cd $JP2DIR
-mkdir merged
+mkdir merge_out
 ls -d */ | parallel -q --jobs $JOBS mergeImages.sh {} "$POLYGON"
 
 find . -name *.xml -type f -delete