enhancement updates to scripts

scripts/skysat_dem_mos.py

@@ -78,7 +78,7 @@ def main():
         print('total dems are {}'.format(len(total_dem_list)))
         out_fn_list = [os.path.join(out_folder,'triplet_{}_mos.tif'.format(stat)) for stat in stats_list]
         print("Mosaicing output total per-pixel nmad, count, nmad and 3 DEMs from 3 stereo combinations in parallel")
-        dem_mos_log = p_map(asp.dem_mosaic,[total_dem_list]*3+[for_aft_list,nadir_aft_list,for_nadir_list],out_fn_list+[os.path.join(out_folder,x) for x in ['for_aft_dem_median_mos.tif', 'nadir_aft_dem_median_mos.tif', 'for_nadir_dem_median_mos.tif']],['None']*6,[None]*6,stats_list+['median']*3,num_cpus=1)
+        dem_mos_log = p_map(asp.dem_mosaic,[total_dem_list]*3+[for_aft_list,nadir_aft_list,for_nadir_list],out_fn_list+[os.path.join(out_folder,x) for x in ['for_aft_dem_median_mos.tif', 'nadir_aft_dem_median_mos.tif', 'for_nadir_dem_median_mos.tif']],['None']*6,[None]*6,stats_list+['median']*3,[None]*6,num_cpus=4)
         out_log_fn = os.path.join(out_folder,'skysat_triplet_dem_mos.log')
         print("Saving triplet DEM mosaic log at {}".format(out_log_fn))
         with open(out_log_fn,'w') as f:
@@ -99,7 +99,7 @@ def main():
         stats_list = ['median','count','nmad']
         print('total dems are {}'.format(len(video_dem_list)))
         out_fn_list = [os.path.join(out_folder,'video_{}_mos.tif'.format(stat)) for stat in stats_list]
-        dem_mos_log = p_map(asp.dem_mosaic,[video_dem_list]*3,out_fn_list,['None']*3,[None]*3,stats_list) 
+        dem_mos_log = p_map(asp.dem_mosaic,[video_dem_list]*3,out_fn_list,['None']*3,[None]*3,stats_list,[None]*3) 
         out_log_fn = os.path.join(out_folder,'skysat_video_dem_mos.log')
         with open(out_log_fn,'w') as f:
             for log in dem_mos_log:

scripts/skysat_orthorectify.py

@@ -82,11 +82,11 @@ def main():
         aft_out_mos = os.path.join(outdir,'aft_map_mos_{}m.tif'.format(tr))
         aft_map_list = sorted(glob.glob(os.path.join(aft_out_dir,'*.tif')))
         print("Preparing forward browse orthomosaic")
-        for_mos_log = asp.dem_mosaic(for_map_list,for_out_mos,tr,tsrs,'first')
+        for_mos_log = asp.dem_mosaic(for_map_list,for_out_mos,tr,tsrs,'first',None)
         print("Preparing nadir browse orthomosaic")
-        nadir_mos_log = asp.dem_mosaic(nadir_map_list, nadir_out_mos, tr, tsrs, 'first')
+        nadir_mos_log = asp.dem_mosaic(nadir_map_list, nadir_out_mos, tr, tsrs, 'first',None)
         print("Preparing aft browse orthomosaic")
-        aft_mos_log = asp.dem_mosaic(aft_map_list, aft_out_mos, tr, tsrs, 'first')
+        aft_mos_log = asp.dem_mosaic(aft_map_list, aft_out_mos, tr, tsrs, 'first',None)
         ## delete temporary files
         if del_opt:
             [shutil.rmtree(x) for x in [for_out_dir,nadir_out_dir,aft_out_dir]]
@@ -150,14 +150,14 @@ def main():
                 median_mosaic = os.path.join(outdir,'{}_{}_{}_median_orthomosaic.tif'.format(for_time,nadir_time,aft_time))
                 wt_avg_mosaic = os.path.join(outdir,'{}_{}_{}_wt_avg_orthomosaic.tif'.format(for_time,nadir_time,aft_time))
                 print("producing finest resolution on top mosaic, per-pixel median and wt_avg mosaic")
-                all_3_view_mos_logs = p_map(asp.dem_mosaic, [res_sorted_list]*3, [res_sorted_mosaic,median_mosaic,wt_avg_mosaic], ['None']*3, [None]*3, ['first','median',None],num_cpus=2)
+                all_3_view_mos_logs = p_map(asp.dem_mosaic, [res_sorted_list]*3, [res_sorted_mosaic,median_mosaic,wt_avg_mosaic], ['None']*3, [None]*3, ['first','median',None],[None]*3,num_cpus=4)
                 res_sorted_log = asp.dem_mosaic(res_sorted_list,res_sorted_mosaic,tr='None',stats='first')
                 print("producing idependent mosaic for different views in parallel")
                 for_mosaic = os.path.join(outdir,'{}_for_first_mosaic.tif'.format(for_time))
                 nadir_mosaic = os.path.join(outdir,'{}_nadir_first_mosaic.tif'.format(nadir_time))
                 aft_mosaic = os.path.join(outdir,'{}_aft_first_mosaic.tif'.format(aft_time))
                 # prepare mosaics in parallel
-                indi_mos_log = p_map(asp.dem_mosaic,[for_img_list,nadir_img_list,aft_img_list], [for_mosaic,nadir_mosaic,aft_mosaic], ['None']*3, [None]*3, ['first']*3,num_cpus=2)
+                indi_mos_log = p_map(asp.dem_mosaic,[for_img_list,nadir_img_list,aft_img_list], [for_mosaic,nadir_mosaic,aft_mosaic], ['None']*3, [None]*3, ['first']*3,[None]*3,num_cpus=4)
                 out_log = os.path.join(outdir,'science_mode_ortho_mos.log')
                 total_mos_log = all_3_view_mos_logs+indi_mos_log
                 print("Saving orthomosaic log at {}".format(out_log))
@@ -173,7 +173,7 @@ def main():
                 print("producing orthomosaic with weighted average statistics")
                 wt_avg_mosaic = os.path.join(outdir,'video_wt_avg_orthomosaic.tif')
                 print("Mosaicing will be done in parallel")
-                all_3_view_mos_logs = p_map(asp.dem_mosaic, [res_sorted_list]*3, [res_sorted_mosaic,median_mosaic,wt_avg_mosaic], ['None']*3, [None]*3, ['first','median',None])
+                all_3_view_mos_logs = p_map(asp.dem_mosaic, [res_sorted_list]*3, [res_sorted_mosaic,median_mosaic,wt_avg_mosaic], ['None']*3, [None]*3, ['first','median',None],[None]*3)
                 out_log = os.path.join(outdir,'science_mode_ortho_mos.log')
                 print("Saving orthomosaic log at {}".format(out_log))
                 with open(out_log,'w') as f:

scripts/skysat_triplet_pipeline.py

@@ -1,13 +1,15 @@
 #! /usr/bin/env python
 import subprocess
 import argparse
-import os,sys,glob
+import os,sys,glob,shutil
 from rpcm import geo
 import numpy as np
 import geopandas as gpd
 from distutils.spawn import find_executable
+from skysat_stereo import misc_geospatial as misc
 from skysat_stereo import asp_utils as asp
 
+
 """
 Script for running the full pipeline based on workflow described in ISPRS 2020 submission
 Need to specify input image folder, input refrence DEM folder
@@ -43,11 +45,14 @@ def main():
     img_list = glob.glob(os.path.join(img_folder,'*.tif'))+glob.glob(os.path.join(img_folder,'*.tiff'))
     if len(img_list)<2:
         print(f"Only {len(img_list)} images detected, exiting")
+        sys.exit()
     if not os.path.exists(coreg_dem):
         print(f"Coreg dem {coreg_dem} could not be located, exiting")
+        sys.exit()
     if not os.path.exists(ortho_dem):
         print(f"Ortho dem {ortho_dem} could not be located, exiting")
-
+        sys.exit()
+
     # structure for output folder
     out_fol = os.path.join(args.outfolder,'proc_out')
     job_name = args.job_name
@@ -102,11 +107,14 @@ def main():
 
     # step 9, final orthorectification
     final_ortho_dir = os.path.join(out_fol,'georegisterd_orthomosaics')
-
-    # step 10, experimental rpc production
+
+    # step 10, plot figure
+    final_figure = os.path.join(out_fol,f"{job_name}_result.jpg")
+
+    # step 11, experimental rpc production
 
     if args.full_workflow == 1:
-         steps2run = np.arange(0,10) # run the entire 9 steps
+         steps2run = np.arange(0,11) # run the entire 9 steps
     else:
         steps2run = np.array(args.partial_workflow_steps).astype(int)
 
@@ -272,6 +280,14 @@ def main():
 		     '-ba_prefix',os.path.join(aligned_cam_dir,'run-run')]
         print("Running final orthomsaic creation")
         asp.run_cmd('skysat_orthorectify.py',ortho_cmd)
+    if 10 in steps2run:
+        # this produces a final plot of orthoimage,DEM, NMAD and countmaps
+        ortho = glob.glob(os.path.join(final_ortho_dir,'*finest_orthomosaic.tif'))[0]
+        count = glob.glob(os.path.join(mos_dem_dir,'*count*.tif'))[0]
+        nmad = glob.glob(os.path.join(mos_dem_dir,'*nmad*.tif'))[0]
+        georegistered_median_dem = glob.glob(os.path.join(alignment_dir,'run-trans_*DEM.tif'))[0]
+        print("plotting final figure")
+        misc.plot_composite_fig(ortho,georegistered_median_dem,count,nmad,outfn=final_figure)
 
 if __name__ == '__main__':
     main()

scripts/skysat_video_pipeline.py

@@ -6,6 +6,7 @@
 import numpy as np
 import geopandas as gpd
 from distutils.spawn import find_executable
+from skysat_stereo import misc_geospatial as misc
 from skysat_stereo import asp_utils as asp
 from skysat_stereo import skysat
 
@@ -100,6 +101,9 @@ def main():
     # step7, final orthorectification
     final_ortho_dir = os.path.join(out_fol,'georegisterd_orthomosaics')
 
+    # step 8, plot figure
+    final_figure = os.path.join(out_fol,f"{job_name}_result.jpg") 
+
     # step 10, experimental rpc production
 
     if args.full_workflow == 1:
@@ -276,5 +280,14 @@ def main():
         print("Running final orthomsaic creation")
         asp.run_cmd('skysat_orthorectify.py',ortho_cmd)
 
+    if 8 in steps2run:
+        # this produces a final plot of orthoimage,DEM, NMAD and countmaps
+        ortho = glob.glob(os.path.join(final_ortho_dir,'*median_orthomosaic.tif'))[0]
+        count = glob.glob(os.path.join(mos_dem_dir,'*count*.tif'))[0]
+        nmad = glob.glob(os.path.join(mos_dem_dir,'*nmad*.tif'))[0]
+        georegistered_median_dem = glob.glob(os.path.join(alignment_dir,'run-trans_*DEM.tif'))[0]
+        print("plotting final figure")
+        misc.plot_composite_fig(ortho,georegistered_median_dem,count,nmad,outfn=final_figure,product='video')
+
 if __name__ == '__main__':
     main()

skysat_stereo/asp_utils.py

@@ -314,8 +314,8 @@ def mapproject(img,outfn,session='rpc',dem='WGS84',tr=None,t_srs='EPSG:4326',cam
         type of input camera model (default: rpc)
     dem: str
         path to input DEM over which images will be draped (default: WGS84, orthorectify just over datum)
-    tr: float/int
-        target resolution of orthorectified output image
+    tr: str
+        target resolution of orthorectified output image (e.g.: '0.9')
     t_srs: str
         target projection of orthorectified output image (default: EPSG:4326)
     cam: str
@@ -357,24 +357,46 @@ def dem_mosaic(img_list,outfn,tr=None,tsrs=None,stats=None,tile_size=None):
         target projection of orthorectified output image (default: EPSG:4326)
     stats: str
         metric to use for mosaicing
+    tile_size: int
+        tile size for distributed mosaicing (if less on memory)
     Returns
     ----------
     out: str
         dem_mosaic log
     """
 
     dem_mosaic_opt = []
-    dem_mosaic_opt.extend(['-o',outfn])
+
     if stats:
         dem_mosaic_opt.extend(['--{}'.format(stats)])
     if (tr is not None) & (ast.literal_eval(tr) is not None):
         dem_mosaic_opt.extend(['--tr', str(tr)])
     if tsrs:
         dem_mosaic_opt.extend(['--t_srs', tsrs])
+    dem_mosaic_args = img_list
     if tile_size:
+        # will first perform tile-wise vertical mosaicing
+        # then blend the result
         dem_mosaic_opt.extend(['--tile-size',str(tile_size)])
-    dem_mosaic_args = img_list
-    out = run_cmd('dem_mosaic',dem_mosaic_args+dem_mosaic_opt)
+        temp_fol = os.path.splitext(outfn)[0]+'_temp'
+        dem_mosaic_opt.extend(['-o',os.path.join(temp_fol,run)])
+        out_tile_op = run_cmd('dem_mosaic',dem_mosaic_args+dem_mosaic_opt)
+        # query all tiles and then do simple mosaic
+        mos_tile_list = sorted(glob.glob(os.path.join(temp_fol+'run-*.tif')))
+        print(f"Found {len(mos_tile_list)}")
+        # now perform simple mosaic
+        dem_mos2_opt = []
+        dem_mos2_opt.extend(['-o',outfn])
+        dem_mos2_args = mos_tile_list
+        out_fn_mos = run_cmd('dem_mosaic',dem_mos2_args+dem_mos2_opt)
+        out = out_tile_op+out_fn_mos
+        print("Deleting tile directory")
+        shutil.rmtree(temp_fol)
+
+    else:
+        # process all at once, no tiling
+        dem_mosaic_opt.extend(['-o',outfn])
+        out = run_cmd('dem_mosaic',dem_mosaic_args+dem_mosaic_opt)
     return out
 
 def get_stereo_opts(session='rpc',threads=4,ba_prefix=None,align='Affineepipolar',xcorr=2,std_mask=0.5,std_kernel=-1,lv=5,corr_kernel=[21,21],rfne_kernel=[35,35],stereo_mode=0,spm=1,cost_mode=2,corr_tile_size=1024,mvs=False):

skysat_stereo/misc_geospatial.py

@@ -1,22 +1,58 @@
 #! /usr/bin/env python
+import matplotlib
+matplotlib.use('Agg')
 
 import numpy as np
 import pandas as pd
 import geopandas as gpd
+from imview import pltlib
+import matplotlib.pyplot as plt
+from pygeotools.lib import iolib,warplib
 
 def shp_merger(shplist):
     """
     merge multiple geopandas shapefiles into 1 multi-row shapefile
     Parameters
-    ----------
+    ------------
     shplist: list
         list of geopandas shapefiles
     Returns
-    ----------
+    ------------
     gpd_merged: geopandas geodataframe
         merged multirow shapefile
     """
     #Taken from here: "https://stackoverflow.com/questions/48874113/concat-multiple-shapefiles-via-geopandas"
     gpd_merged = pd.concat([shp for shp in shplist]).pipe(gpd.GeoDataFrame)
     return gpd_merged
 
+def plot_composite_fig(ortho,dem,count,nmad,outfn,product='triplet'):
+    """
+    Plot the gallery figure for final DEM products
+    Parameters
+    ------------
+    ortho: str
+        path to orthoimage
+    dem: str
+        path to dem
+    count: str
+        path to count map
+    nmad: str
+        path to NMAD
+    outfn: str
+        path to save output figure
+    ortho: str
+        product to plot (triplet/video)
+    """
+    if product == 'triplet':
+        figsize=(10,8)
+    else:
+        figsize=(10,3)
+    f,ax = plt.subplots(1,4,figsize=figsize)
+    ds_list = warplib.memwarp_multi_fn([ortho,dem,count,nmad],res='max')
+    ortho,dem,count,nmad = [iolib.ds_getma(x) for x in ds_list]
+    pltlib.iv(ortho,ax=ax[0],cmap='gray',scalebar=True,cbar=False,ds=ds_list[0],skinny=False)
+    pltlib.iv(dem,ax=ax[1],hillshade=True,scalebar=False,ds=ds_list[1],label='Elevation (m WGS84)',skinny=False)
+    pltlib.iv(count,ax=ax[2],cmap='YlOrRd',label='DEM count',skinny=False)
+    pltlib.iv(nmad,ax=ax[3],cmap='inferno',clim=(0,10),label='Elevation NMAD (m)',skinny=False)
+    plt.tight_layout()
+    f.savefig(outfn,dpi=300,bbox_inches='tight',pad_inches=0.1)