diff --git a/Common/geometry.py b/Common/geometry.py
index 6d7397e9f97bd59ceb38ec0dc5de8b6abacf0f6d..235b30ed8b5ef60d0855417fac3c6943b807accd 100644
--- a/Common/geometry.py
+++ b/Common/geometry.py
@@ -1,3 +1,4 @@
+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
+
diff --git a/Common/imagetools.py b/Common/imagetools.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a0900fd2ed4ac23d7c0b0f9a24d8f1c1e05fa79
--- /dev/null
+++ b/Common/imagetools.py
@@ -0,0 +1,18 @@
+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])
diff --git a/Workflows/operations.py b/Workflows/operations.py
index a80026c680a3d2d9620876e50f9d95c59332055e..a861c2e0f8dd5720a7391e5ed0f13683377c9026 100644
--- a/Workflows/operations.py
+++ b/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
diff --git a/moringa.py b/moringa.py
index c2a8b2189424ce31de47108da95c98e027bec439..771c3b52fed8f068645beb114d34a84149da8adc 100755
--- a/moringa.py
+++ b/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