UDPATE: incidence locale

python/SoilMoisturePipeline.py

@@ -1030,7 +1030,13 @@ if __name__ == "__main__":
                 sarth = os.path.join(sard , "incidenceAngleFromEllipsoid.img")
             else:
                 sarvv = sard
-                sarth = sard[:-6]+"THETA.TIF"
+                sarbase = os.path.dirname(sarvv)
+                if os.path.exists(os.path.join(sarbase,"auxfiles")) :
+                    auxdir = os.path.join(sarbase,"auxfiles")
+                    thfname = os.path.basename(sard).replace("_VV.TIF","_THETALOC.TIF")
+                    sarth = os.path.join(auxdir,thfname)
+                else:
+                    sarth = sard[:-6]+"THETA.TIF"
             
             if args.outformat == "raster":
                 output = os.path.join(args.outdir, "MV_"+sarsplit[0]+"_"+args.zone+"_"+sarsplit[4]+".TIF") 

python/SoilMoistureS1Calibration.py

@@ -21,7 +21,6 @@ namespaces = {
     'gml': "http://www.opengis.net/gml"
     }
 overwrite=False
-slope=True
 
 def search_files(directory='.', resolution='S1', extension='SAFE', fictype='d'):
     images=[]
@@ -181,7 +180,7 @@ def prepare_srtm_hgt(sarsafe):
     print("Ok for SRTM tiles: "+str(tilesname))
     return tiftiles
 
-def create_zone_dem_and_slope(ref,zone,outdir):
+def create_zone_dem_aspect_and_slope(ref,zone,outdir):
     
     outdem = os.path.join(outdir,"SRTM_"+zone+".TIF")
     if os.path.exists(outdem) and not overwrite :
@@ -210,17 +209,22 @@ def create_zone_dem_and_slope(ref,zone,outdir):
     supapp.SetParameterString("out", outdem)
     supapp.ExecuteAndWriteOutput()
     
-    #Compute slope
-    if slope :
-        print("Generating slope file...")
-        outslope = os.path.join(outdir,"SLOPE_"+zone+".TIF")
-        cmd=['gdaldem', 'slope', '-p', outdem, outslope]
-        
-        process_command(cmd)
+    #Compute slope & aspect
+    print("Generating slope file...")
+    outslope = os.path.join(outdir,"SLOPE_"+zone+".TIF")
+    cmd=['gdaldem', 'slope', '-p', outdem, outslope]
+    
+    process_command(cmd)
+    
+    print("Generating aspect file...")
+    outaspect = os.path.join(outdir,"ASPECT_"+zone+".TIF")
+    cmd=['gdaldem', 'aspect', outdem, outaspect]
+    
+    process_command(cmd)
     
     return outdem
         
-def create_incidenceAngle_raster(annofiles,dem, img):
+def create_incidenceAngle_from_surface_raster(annofiles,dem, img):
     
     outdir = os.path.dirname(img)
     tempdir = os.path.join(os.path.dirname(img),"auxfiles")
@@ -233,8 +237,8 @@ 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 os.path.exists(outinc) and not overwrite :
+        return outinc
     
     if not os.path.exists(tempdir):
         os.mkdir(tempdir)
@@ -284,11 +288,24 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     process_command(cmd)
     
     
+    return outinc
+
+def create_incidenceAngle_from_ellipsoid_raster(vvortho,dem,thetasurf):
+    
+    outdir = os.path.dirname(thetasurf)
+    
+    outthetaelli = thetasurf.replace("_THETASURF","_THETAELLI")  
+    
+    
+    if os.path.exists(outthetaelli) and not overwrite :
+        return outthetaelli
+    
+    
     app41 = otbApplication.Registry.CreateApplication("Superimpose")
     # The following lines set all the application parameters:
-    app41.SetParameterString("inm", outinc)
-    app41.SetParameterString("inr", img)
-    app41.SetParameterString("interpolator","nn")
+    app41.SetParameterString("inm", thetasurf)
+    app41.SetParameterString("inr", vvortho)
+    app41.SetParameterString("interpolator","bco")
     app41.SetParameterString("out", "temp41.tif")
     
     print("Launching... Resampling")
@@ -299,8 +316,8 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     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("inr", vvortho)
+    app42.SetParameterString("interpolator","bco")
     app42.SetParameterString("out", "temp41.tif")
     
     print("Launching... Resampling")
@@ -309,7 +326,7 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     print("End of Resampling \n")
     
     appS = otbApplication.Registry.CreateApplication("BandMath")
-    appS.SetParameterStringList("il", [img])
+    appS.SetParameterStringList("il", [vvortho])
     appS.AddImageToParameterInputImageList("il", app42.GetParameterOutputImage("out"))
     appS.AddImageToParameterInputImageList("il", app41.GetParameterOutputImage("out"))
     # Define Input im2: Band Red (B4)
@@ -320,11 +337,101 @@ def create_incidenceAngle_raster(annofiles,dem, img):
     
     appM = otbApplication.Registry.CreateApplication("ManageNoData")
     appM.SetParameterInputImage("in", appS.GetParameterOutputImage("out"))
-    appM.SetParameterString("out", outincfromdem+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
+    appM.SetParameterString("out", outthetaelli+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
     appM.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
     appM.SetParameterString("mode", "changevalue")
     appM.ExecuteAndWriteOutput()
     
+    return outthetaelli
+
+def create_local_incidenceAngle_raster(annofiles, img, dem, thetasurf):
+    
+    outdir = os.path.dirname(img)
+    
+    inslope = dem.replace("SRTM_", "SLOPE_")
+    inaspect = dem.replace("SRTM_", "ASPECT_")
+    #todo: mettre un dosiier temp
+    outthetaloc = thetasurf.replace("_THETASURF","_THETALOC")  
+    
+    if os.path.exists(outthetaloc) and not overwrite :
+        return outthetaloc
+    
+    
+    plateformheading=[]
+    for annofile in annofiles:
+    
+        tree = ET.parse(annofile)
+        root = tree.getroot()
+        
+        plateformheading.append(float(root.find(".//platformHeading").text))
+        
+    if len(plateformheading) == 0:
+        print("Error: Can't find plateformHeading in annotation files")
+        exit()
+    elif len(plateformheading) == 1:
+        azimuthsat = (plateformheading[0]+90) % 360
+    else:
+        azimuthsat = (math.fsum(plateformheading)/len(plateformheading)+90) % 360
+        
+
+    
+    app41 = otbApplication.Registry.CreateApplication("Superimpose")
+    # The following lines set all the application parameters:
+    app41.SetParameterString("inm", thetasurf)
+    app41.SetParameterString("inr", img)
+    app41.SetParameterString("interpolator","bco")
+    app41.SetParameterString("out", "temp41.tif")
+    
+    print("Launching... Resampling")
+    # The following line execute the application
+    app41.Execute()
+    print("End of Resampling \n")
+    
+    app43 = otbApplication.Registry.CreateApplication("Superimpose")
+    # The following lines set all the application parameters:
+    app43.SetParameterString("inm", inslope)
+    app43.SetParameterString("inr", img)
+    app43.SetParameterString("interpolator","bco")
+    app43.SetParameterString("out", "temp43.tif")
+    
+    print("Launching... Resampling")
+    # The following line execute the application
+    app43.Execute()
+    print("End of Resampling \n")
+    
+    app44 = otbApplication.Registry.CreateApplication("Superimpose")
+    # The following lines set all the application parameters:
+    app44.SetParameterString("inm", inaspect)
+    app44.SetParameterString("inr", img)
+    app44.SetParameterString("interpolator","bco")
+    app44.SetParameterString("out", "temp44.tif")
+    
+    print("Launching... Resampling")
+    # The following line execute the application
+    app44.Execute()
+    print("End of Resampling \n")
+    
+    appS = otbApplication.Registry.CreateApplication("BandMath")
+    appS.SetParameterStringList("il", [img])  #im1b1 -> vv
+    appS.AddImageToParameterInputImageList("il", app41.GetParameterOutputImage("out")) #im2b1 -> incidence surf
+    appS.AddImageToParameterInputImageList("il", app43.GetParameterOutputImage("out")) #im3b1 -> slope
+    appS.AddImageToParameterInputImageList("il", app44.GetParameterOutputImage("out")) #im4b1 -> aspect
+    # Define Input im2: Band Red (B4)
+    appS.SetParameterString("out", "tempS.tif")
+    appS.SetParameterString("exp", "im1b1 == 0?0:acos(cos(im3b1*"+str(math.pi/400)+")*cos(im2b1*"+str(math.pi/180)+")-sin(im3b1*"+str(math.pi/400)+")*sin(im2b1*"+str(math.pi/180)+")*cos(("+str(azimuthsat)+"-im4b1)*"+str(math.pi/180)+"))*"+str(180/math.pi)+"")
+    
+    appS.Execute()
+    
+    appM = otbApplication.Registry.CreateApplication("ManageNoData")
+    appM.SetParameterInputImage("in", appS.GetParameterOutputImage("out"))
+    appM.SetParameterString("out", outthetaloc+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
+    appM.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
+    appM.SetParameterString("mode", "changevalue")
+    appM.ExecuteAndWriteOutput()
+    
+    return outthetaloc
+
+    
     
 def cal_ortho_sub_mosa_pipeline(sarfiles, ref, outfile):
     
@@ -392,29 +499,27 @@ if __name__ == "__main__":
     # Make parser object
     parser = argparse.ArgumentParser(description=
         """
-        S1 Calibration, orhorectification, dem, slope and incidence angle processing
+        S1 Calibration, orhorectification, dem, slope, aspect and incidence angle processing
         """)
     
     parser.add_argument('-indir', action='store', required=True, help='Directory containing Sentinel1 SAFE dir')
     parser.add_argument('-inref', action='store', required=True, help='Image reference, generally corresponding to a S2 tile')
     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('--no-slope', dest='slope', action='store_false',required=False, help="[Optional] Don't create slope from dem (default True), MUST have ref image in projected coordinates")
+    parser.add_argument('-theta', choices=['loc', 'elli'],  default='loc', required=False, help='[Optional]Create local incidence raster or incidence raster from ellipsoid, default loc')
     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()
     
-    slope = args.slope
     overwrite = args.overwrite
     
-    if slope:
-        raster = gdal.Open(args.inref)
-        proj = osr.SpatialReference(wkt=raster.GetProjection())
-        if proj.IsGeographic() :
-            print("Error: Reference image must be in projected coordinates in order to create slope")
-            exit()
+    raster = gdal.Open(args.inref)
+    proj = osr.SpatialReference(wkt=raster.GetProjection())
+    if proj.IsGeographic() :
+        print("Error: Reference image must be in projected coordinates in order to create slope")
+        exit()
     
     if not os.path.exists(args.outdir):
         os.mkdir(args.outdir)
@@ -427,7 +532,7 @@ if __name__ == "__main__":
     indir=search_files(args.indir, 'S1', 'SAFE', 'd')
     
     
-    zonedemfile = create_zone_dem_and_slope(args.inref,args.zone,demdir)
+    zonedemfile = create_zone_dem_aspect_and_slope(args.inref,args.zone,demdir)
     
     dejafait=[]
     for file in indir:
@@ -510,8 +615,10 @@ if __name__ == "__main__":
         else :
             print("Error on index :"+l)
     
-        
-        
-        create_incidenceAngle_raster(annofiles,zonedemfile,outfile)
+        thetasurf = create_incidenceAngle_from_surface_raster(annofiles,zonedemfile, outfile)
+        if args.theta == "elli":
+            create_incidenceAngle_from_ellipsoid_raster(outfile,zonedemfile,thetasurf)
+        elif args.theta == "loc":
+            create_local_incidenceAngle_raster(annofiles, outfile, zonedemfile, thetasurf)