ENH: added tools for rigid shift and basemap creation.

Common/geometry.py

+import os
 import sys
 
 import numpy as np
@@ -10,6 +11,8 @@ from typing import List
 from tqdm import tqdm
 from datetime import *
 import multiprocessing as mp
+import pickle
+import rasterio
 
 def compute_overlap_matrix(_lst: List[otb.Application], out_param='out'):
     masks = []
@@ -388,3 +391,19 @@ def get_patch_centers(msk, patch_size, n_patches, min_coverage=0, min_cov_extent
         return np.asarray(out_coords), coverage, cov_extent
     else:
         return None, coverage, cov_extent
+
+def apply_shift_using_points(img, pts_file, separator=','):
+    with open(pts_file, 'r') as pts_f:
+        pts = pts_f.read().splitlines()
+    pts = [x.split(separator) for x in pts if x]
+    pts = [[float(x) for x in y] for y in pts]
+    pts = np.array(pts)
+    mns = np.array([pts[:,2]-pts[:,0],pts[:,3]-pts[:,1]])
+    mns = np.mean(mns, axis=1)
+    with rasterio.open(img, 'r+') as ds:
+        tr = ds.transform
+        with open(os.path.splitext(img)[0]+'_transform_backup.pkl', 'wb') as f:
+            pickle.dump(tr, f)
+        ds.transform = rasterio.transform.Affine(tr[0], tr[1], tr[2] + mns[0], tr[3], tr[4], tr[5] + mns[1])
+    return mns
+    

Common/imagetools.py

+import os
+import otbApplication as otb
+import rasterio
+
+def vhr_to_basemap(img,out,quantiles=[2,2], bands=[1,2,3], quality=75):
+    dc = otb.Registry.CreateApplication('DynamicConvert')
+    dc.SetParameterString('in',img)
+    dc.SetParameterInt('quantile.low', quantiles[0])
+    dc.SetParameterInt('quantile.high', quantiles[1])
+    dc.SetParameterString('channels', 'rgb')
+    dc.SetParameterInt('channels.rgb.red', bands[0])
+    dc.SetParameterInt('channels.rgb.green', bands[1])
+    dc.SetParameterInt('channels.rgb.blue', bands[2])
+    dc.SetParameterString('out', out + '?&gdal:co:compress=jpeg&gdal:co:jpeg_quality={}'.format(quality))
+    dc.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint8)
+    dc.ExecuteAndWriteOutput()
+    with rasterio.open(out,'r+') as ds:
+        ds.build_overviews(factors=[2,4,8,16,32,64])

Workflows/operations.py

@@ -3,6 +3,8 @@ import OBIA.segmentation
 import VHR.vhrbase
 from TimeSeries import s2theia, s2planetary, s1base, s1planetary, planet_mosaics, landsat_planetary
 from Common import demtools
+from Common.geometry import apply_shift_using_points
+from Common.imagetools import vhr_to_basemap
 
 def run_segmentation(img, threshold, cw, sw , out_seg,
                      n_first_iter, margin, roi, n_proc, memory=None,
@@ -124,3 +126,11 @@ def fetch(imagery, shp, out_fld, dt=None, auth=None, only_tiles=None):
     elif imagery == 'landsatplanetary':
         landsat_planetary.fetch(shp, dt, out_fld, auth)
     return
+
+def apply_shift(in_img, in_points):
+    apply_shift_using_points(in_img, in_points)
+    return
+
+def create_basemap(img, out, quantiles=[2,2], bands=[1,2,3], quality=75):
+    vhr_to_basemap(img, out, quantiles, bands, quality)
+    return
\ No newline at end of file

moringa.py

@@ -90,6 +90,23 @@ def main(args):
     dategen.add_argument("end_date", type=str, help="End date in YYYY-MM-DD format.")
     dategen.add_argument("step", type=int, help="Number of days between two dates.")
 
+    applyshift = subpar.add_parser("apply_shift", help="Apply a rigid geotransform shift using a set of tie points.",
+                                   formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    applyshift.add_argument("input_image", type=str, help="Path to the input image. Shift is applied in-place.")
+    applyshift.add_argument("input_points", type=str, help="Path to the tie points file (x_img,y_img,x_ref,y_ref)")
+
+    basemap = subpar.add_parser("create_basemap", help="Generate a browsing optimized version of an image.",
+                                   formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    basemap.add_argument("input_image", type=str, help="Path to the input image.")
+    basemap.add_argument("output_image", type=str, help="Path to the tie points file (x_img,y_img,x_ref,y_ref)")
+    basemap.add_argument("--low_quantile", type=int, default=2, help="Low quantile for dynamic conversion.")
+    basemap.add_argument("--high_quantile", type=int, default=2, help="High quantile for dynamic conversion.")
+    basemap.add_argument("--red", type=int, default=1, help="Band for red channel.")
+    basemap.add_argument("--green", type=int, default=2, help="Band for green channel.")
+    basemap.add_argument("--blue", type=int, default=3, help="Band for blue channel.")
+    basemap.add_argument("--quality", type=int, default=75, help="Quality (0-100) for jpeg conversion.")
+
+
     if len(args) == 1:
         parser.print_help()
         sys.exit(0)
@@ -127,6 +144,13 @@ def main(args):
 
     if arg.cmd == "gen_date_file":
         periodic_date_generator(arg.output_file, arg.start_date, arg.end_date, arg.step)
+    
+    if arg.cmd == "apply_shift":
+        apply_shift(arg.input_image, arg.input_points)
+
+    if arg.cmd == "create_basemap":
+        create_basemap(arg.input_image, arg.output_image, quantiles=[arg.low_quantile, arg.high_quantile],
+                       bands=[arg.red,arg.green,arg.blue], quality=arg.quality)
 
     return 0