UPDATE: compress theta + overwrite option

python/SoilMoistureS1Calibration.py

@@ -20,7 +20,7 @@ globconfig = {
 namespaces = {
     'gml': "http://www.opengis.net/gml"
     }
-
+overwrite=False
 
 def search_files(directory='.', resolution='S1', extension='SAFE', fictype='d'):
     images=[]
@@ -182,6 +182,10 @@ def prepare_srtm_hgt(sarsafe):
 
 def create_zone_dem_and_slope(ref,zone,outdir):
     
+    outdem = os.path.join(outdir,"SRTM_"+zone+".TIF")
+    if os.path.exists(outdem) and not overwrite :
+        return outdem
+        
     #get tiles
     pointLL, pointUR, cols, rows = get_img_extend(ref)
     
@@ -193,7 +197,7 @@ def create_zone_dem_and_slope(ref,zone,outdir):
         tilefiles.append(os.path.join(tifdir,tile+".tif"))
     
     print("Generating DEM file on ref:"+ref)
-    outdem = os.path.join(outdir,"SRTM_"+zone+".TIF")
+    
     mosapp = otbApplication.Registry.CreateApplication("Mosaic")
     mosparams = {"il":tilefiles, "out":"mosatemp.tif"}
     mosapp.SetParameters(mosparams)
@@ -227,47 +231,50 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     csvfile = os.path.join(tempdir,"_".join(imgsp[:-1])+"_THETA.csv")
     layername = "_".join(imgsp[:-1])+"_THETA"
     
+    if os.path.exists(outincfromdem) and not overwrite :
+        return outincfromdem
+    
     if not os.path.exists(tempdir):
         os.mkdir(tempdir)
     
     pointLL, pointUR, cols, rows = get_img_extend(img)
     
+    if not os.path.exists(csvfile):
+        for annofile in annofiles:
     
-    for annofile in annofiles:
-
-        tree = ET.parse(annofile)
-        root = tree.getroot()
-        
-        geoloclist = root.find(".//geolocationGridPointList")
-        
-        if geoloclist == None :
-            print("Can't find geolocationGridPointList tag in xml file.")
-            exit()
-        
-        header = True
-        if os.path.exists(csvfile):
-            header = False
-        if not os.path.exists(vrtfile):
-            with open(vrtfile,'w') as vrt:
-                vrt.write("<OGRVRTDataSource>\n")
-                vrt.write('    <OGRVRTLayer name="'+layername+'">\n')
-                vrt.write('        <SrcDataSource>'+csvfile+'</SrcDataSource>\n')
-                vrt.write("        <GeometryType>wkbPoint</GeometryType>\n")
-                vrt.write("        <LayerSRS>WGS84</LayerSRS>\n")
-                vrt.write('        <GeometryField encoding="PointFromColumns" x="Easting" y="Northing" z="incidence"/>\n')
-                vrt.write("    </OGRVRTLayer>\n")
-                vrt.write("</OGRVRTDataSource>\n")
+            tree = ET.parse(annofile)
+            root = tree.getroot()
+            
+            geoloclist = root.find(".//geolocationGridPointList")
             
-        with open(csvfile,'a') as csvf:
-            if header:
-                csvf.write("Easting,Northing,incidence\n")
-            for child in geoloclist:
-                lat = child.find(".//latitude").text
-                lon = child.find(".//longitude").text
-                alti = float(child.find(".//height").text)
-                incdem = float(child.find(".//incidenceAngle").text)
-                inc = incdem - alti/693000.
-                csvf.write(lon+","+lat+","+str(incdem)+"\n")
+            if geoloclist == None :
+                print("Can't find geolocationGridPointList tag in xml file.")
+                exit()
+            
+            header = True
+            if os.path.exists(csvfile):
+                header = False
+            if not os.path.exists(vrtfile):
+                with open(vrtfile,'w') as vrt:
+                    vrt.write("<OGRVRTDataSource>\n")
+                    vrt.write('    <OGRVRTLayer name="'+layername+'">\n')
+                    vrt.write('        <SrcDataSource>'+csvfile+'</SrcDataSource>\n')
+                    vrt.write("        <GeometryType>wkbPoint</GeometryType>\n")
+                    vrt.write("        <LayerSRS>WGS84</LayerSRS>\n")
+                    vrt.write('        <GeometryField encoding="PointFromColumns" x="Easting" y="Northing" z="incidence"/>\n')
+                    vrt.write("    </OGRVRTLayer>\n")
+                    vrt.write("</OGRVRTDataSource>\n")
+                
+            with open(csvfile,'a') as csvf:
+                if header:
+                    csvf.write("Easting,Northing,incidence\n")
+                for child in geoloclist:
+                    lat = child.find(".//latitude").text
+                    lon = child.find(".//longitude").text
+                    alti = float(child.find(".//height").text)
+                    incdem = float(child.find(".//incidenceAngle").text)
+                    inc = incdem - alti/693000.
+                    csvf.write(lon+","+lat+","+str(incdem)+"\n")
      
     
     cmd = ["gdal_grid","-a","linear","-txe", str(pointLL.GetX()), str(pointUR.GetX()), "-tye", str(pointLL.GetY()), str(pointUR.GetY()), "-outsize", str(cols), str(rows), "-of", "GTiff", "-ot", "Float32", "-l", layername, vrtfile, outinc]
@@ -278,7 +285,7 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     app41 = otbApplication.Registry.CreateApplication("Superimpose")
     # The following lines set all the application parameters:
     app41.SetParameterString("inm", outinc)
-    app41.SetParameterString("inr", dem)
+    app41.SetParameterString("inr", img)
     app41.SetParameterString("interpolator","nn")
     app41.SetParameterString("out", "temp41.tif")
     
@@ -287,14 +294,34 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     app41.Execute()
     print("End of Resampling \n")
     
+    app42 = otbApplication.Registry.CreateApplication("Superimpose")
+    # The following lines set all the application parameters:
+    app42.SetParameterString("inm", dem)
+    app42.SetParameterString("inr", img)
+    app42.SetParameterString("interpolator","nn")
+    app42.SetParameterString("out", "temp41.tif")
+    
+    print("Launching... Resampling")
+    # The following line execute the application
+    app42.Execute()
+    print("End of Resampling \n")
+    
     appS = otbApplication.Registry.CreateApplication("BandMath")
-    appS.SetParameterStringList("il", [dem])
+    appS.SetParameterStringList("il", [img])
+    appS.AddImageToParameterInputImageList("il", app42.GetParameterOutputImage("out"))
     appS.AddImageToParameterInputImageList("il", app41.GetParameterOutputImage("out"))
     # Define Input im2: Band Red (B4)
-    appS.SetParameterString("out", outincfromdem)
-    appS.SetParameterString("exp", "im2b1 - im1b1/693000" )
+    appS.SetParameterString("out", "tempS.tif")
+    appS.SetParameterString("exp", "im1b1 == 0?0:im3b1 - im2b1/693000" )
     
-    appS.ExecuteAndWriteOutput()
+    appS.Execute()
+    
+    appM = otbApplication.Registry.CreateApplication("ManageNoData")
+    appM.SetParameterInputImage("in", appS.GetParameterOutputImage("out"))
+    appM.SetParameterString("out", outincfromdem+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
+    appM.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
+    appM.SetParameterString("mode", "changevalue")
+    appM.ExecuteAndWriteOutput()
     
     
 def cal_ortho_sub_mosa_pipeline(sarfiles, ref, outfile):
@@ -370,12 +397,14 @@ if __name__ == "__main__":
     parser.add_argument('-inref', action='store', required=True, help='Image references Sentinel-2 Tile, example any S2 NDVI')
     parser.add_argument('-zone', action='store', required=True, help='Geographic zone reference for output files naming')
     parser.add_argument('-outdir', action='store', required=True, help='Output directory')
-#    parser.add_argument('-user', action='store', required=True, help='SciHub Copernicus user name')
-#    parser.add_argument('-psswd', action='store', required=True, help='SciHub Copernicus user password')
+    parser.add_argument('-overwrite', dest='overwrite', action='store_true',required=False, help='[Optional] Overwrite already existing files (default False)')
+    parser.set_defaults(overwrite=False)
     
     
     args=parser.parse_args()
     
+    overwrite = args.overwrite
+    
     if not os.path.exists(args.outdir):
         os.mkdir(args.outdir)
         
@@ -386,6 +415,7 @@ if __name__ == "__main__":
     indir=[]
     indir=search_files(args.indir, 'S1', 'SAFE', 'd')
     
+    
     zonedemfile = create_zone_dem_and_slope(args.inref,args.zone,demdir)
     
     dejafait=[]
@@ -420,7 +450,7 @@ if __name__ == "__main__":
             dejafait.append(indir[l[0]])
             dejafait.append(indir[l[1]])
             
-            if not os.path.exists(outfile) :        
+            if not os.path.exists(outfile) or overwrite :        
                 cal_ortho_sub_mosa_pipeline(sarfiles, args.inref, outfile)
             
             
@@ -450,7 +480,7 @@ if __name__ == "__main__":
             
             outfile = outfilebase +"_VV.TIF"
             
-            if not os.path.exists(outfile) :        
+            if not os.path.exists(outfile) or overwrite :         
                 cal_ortho_sub_mosa_pipeline(in1vv, args.inref, outfile)