diff --git a/Common/geometry.py b/Common/geometry.py
index 235b30ed8b5ef60d0855417fac3c6943b807accd..365ab8ac0b625d39434f2109aee34cc090a76a49 100644
--- a/Common/geometry.py
+++ b/Common/geometry.py
@@ -13,6 +13,7 @@ from datetime import *
import multiprocessing as mp
import pickle
import rasterio
+from scipy.stats import gaussian_kde
def compute_overlap_matrix(_lst: List[otb.Application], out_param='out'):
masks = []
@@ -49,11 +50,24 @@ def local_shift(src, tgt):
return shifts
+def compute_mode_2d(arr):
+ sm = gaussian_kde(arr.transpose())
+ mx, Mx = min(arr[:,0]), max(arr[:,0])
+ my, My = min(arr[:,1]), max(arr[:,1])
+ X, Y = np.mgrid[mx:Mx:0.1, my:My:0.1]
+ positions = np.vstack([X.ravel(), Y.ravel()])
+ Z = np.reshape(sm(positions).T, X.shape)
+ idx = np.unravel_index(Z.argmax(), Z.shape)
+ out = X[idx[0],0], Y[0,idx[1]]
+ return np.array(out)
+
+
+
def compute_displacement(_src: otb.Application, _tgt: otb.Application,
src_band=2, tgt_band=2,
out_param_src='out', out_param_tgt='out',
geobin_size=32, geobin_spacing=256, margin=32,
- filter=5, n_proc=6):
+ filter=5, n_proc=6, use_kde=False):
sz = _tgt.GetImageSize(out_param_tgt)
W,H = sz[0], sz[1]
@@ -102,7 +116,10 @@ def compute_displacement(_src: otb.Application, _tgt: otb.Application,
shifts = shifts[nrm < filter]
if len(shifts) > 0:
- return np.mean(np.array(shifts), axis=0), len(shifts)
+ if not use_kde:
+ return np.mean(np.array(shifts), axis=0), len(shifts)
+ else:
+ return compute_mode_2d(shifts), len(shifts)
else:
return None, 0
diff --git a/Common/otb_numpy_proc.py b/Common/otb_numpy_proc.py
index b784ce42a0a0b891495e33b2ccc824af43e18114..94adf8feb42f58c4e8bebbd48aafa20f6166bec3 100644
--- a/Common/otb_numpy_proc.py
+++ b/Common/otb_numpy_proc.py
@@ -5,6 +5,8 @@ import os, psutil
from functools import partial
import string, random
from osgeo.gdal import BuildVRT
+from osgeo.osr import SpatialReference
+from pyproj import Transformer as T
import otbApplication as otb
def randomword(length):
@@ -205,6 +207,26 @@ def to_otb_pipeline(obj):
sys.exit("Impossible to convert object ot OTB pipeline")
return ppl
+def ppl_output_extent(img, out_srs=None, param='out'):
+ og = list(img.GetImageOrigin(param))
+ sp = list(img.GetImageSpacing(param))
+ sz = list(img.GetImageSize(param))
+ xt = [og[0], og[0]+sp[0]*sz[0], og[1], og[1]+sp[1]*sz[1]]
+ if xt[0] > xt[1]:
+ xt[0], xt[1] = xt[1], xt[0]
+ if xt[2] > xt[3]:
+ xt[2], xt[3] = xt[3], xt[2]
+ if out_srs==None:
+ return xt
+ else:
+ i_srs = SpatialReference()
+ i_srs.ImportFromWkt(img.GetImageProjection(param))
+ tr = T.from_crs(i_srs.GetAuthorityCode(None), out_srs, always_xy=True)
+ lon_min, lat_min = tr.transform(xt[0], xt[2])
+ lon_max, lat_max = tr.transform(xt[1], xt[3])
+ bbox = [lon_min, lon_max, lat_min, lat_max]
+ return bbox
+
def extract_roi_as_numpy(img: otb.Application, startx=None, starty=None, sizex=None, sizey=None, bands=None, out_param='out'):
ers = otb.Registry.CreateApplication('ExtractROI')
ers.SetParameterInputImage('in', img.GetParameterOutputImage(out_param))
diff --git a/VHR/googlesat.py b/VHR/googlesat.py
index 9aecd1b9548e6ee5c620da3a9135240232fc6947..799261210b85721cd389a1d00c209134f40fe919 100644
--- a/VHR/googlesat.py
+++ b/VHR/googlesat.py
@@ -4,7 +4,7 @@ import os
from tqdm import tqdm
from osgeo import gdal
import otbApplication as otb
-from Common.otb_numpy_proc import to_otb_pipeline
+from Common.otb_numpy_proc import to_otb_pipeline, ppl_output_extent
class GoogleSatScene:
@@ -31,6 +31,9 @@ class GoogleSatScene:
"{y}", str(y)).replace(
"{z}", str(z))
path = f'{temp_dir}/{x}_{y}_{z}.png'
+ opener = urllib.request.build_opener()
+ opener.addheaders = [('User-agent', 'Mozilla/5.0')]
+ urllib.request.install_opener(opener)
urllib.request.urlretrieve(url, path)
return(path)
@@ -67,7 +70,15 @@ class GoogleSatScene:
outputBounds=bounds)
return filename + '.tif'
- def __init__(self, bbox, out, zoom=18):
+ def __init__(self, roi, out, zoom=18):
+
+ if isinstance(roi, list):
+ bbox = roi
+ elif isinstance(roi, otb.Application):
+ #if ROI given with otbApplication object, image parameter must be 'out'
+ bbox = ppl_output_extent(roi, out_srs=4326)
+
+ print(bbox)
fld = os.path.dirname(out)
s = os.path.basename(out)
@@ -85,14 +96,28 @@ class GoogleSatScene:
tile_lst.append(self.georeference_raster_tile(x, y, zoom, png_path))
os.remove(png_path)
pbar.update()
+
+ pbar.close()
mos = otb.Registry.CreateApplication('Mosaic')
mos.SetParameterStringList('il', tile_lst)
- mos.SetParameterString('out', out)
- mos.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint8)
- mos.ExecuteAndWriteOutput()
+
+
+ if isinstance(roi, list):
+ mos.SetParameterString('out', out)
+ mos.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint8)
+ mos.ExecuteAndWriteOutput()
+ elif isinstance(roi, otb.Application):
+ mos.Execute()
+ si = otb.Registry.CreateApplication('Superimpose')
+ si.SetParameterInputImage('inr', roi.GetParameterOutputImage('out'))
+ si.SetParameterInputImage('inm', mos.GetParameterOutputImage('out'))
+ si.SetParameterString('out', out)
+ si.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint8)
+ si.ExecuteAndWriteOutput()
+
ds = gdal.Open(out,1)
[ds.GetRasterBand(i).DeleteNoDataValue() for i in range(1, ds.RasterCount+1)]
-
- [os.remove(x) for x in tile_lst]
self.scene = to_otb_pipeline(out)
+ [os.remove(x) for x in tile_lst]
+ return
\ No newline at end of file
diff --git a/VHR/vhrbase.py b/VHR/vhrbase.py
index 89f58ed20f02f3720a028fa9ccc142bbd6c7657e..88efccdc2828826b9a2feeec7dbcc9094eb0289b 100644
--- a/VHR/vhrbase.py
+++ b/VHR/vhrbase.py
@@ -7,6 +7,8 @@ import re
from osgeo import ogr, gdal
from Common.geometry import compute_displacement
from math import ceil, floor, sqrt
+import rasterio
+import pickle
class SPOT67RasterPipeline:
# BEGIN SPOT6/7 VHR PROTOTYPE
@@ -272,4 +274,37 @@ class PleiadesOrthoAutoMosaicPipeline:
return fn
else:
print('Preprocess scenes first.')
- return None
\ No newline at end of file
+ return None
+
+class PreprocessedVHR:
+ def __init__(self, path, red_channel=2):
+ self.path = path
+ self.scene = to_otb_pipeline(path)
+ self.spacing = self.scene.GetImageSpacing('out')
+ self.red_channel = red_channel
+
+ def shift_from_reference(self, ref, red_channel=2,
+ geobin_size=32, geobin_spacing=32, margin=32, filter=20):
+ rfp = otb.Registry.CreateApplication('Superimpose')
+ rfp.SetParameterInputImage('inr', self.scene.GetParameterOutputImage('out'))
+ rfp.SetParameterString('inm', ref)
+ rfp.Execute()
+ d, n = compute_displacement(self.scene, rfp, self.red_channel, red_channel,
+ geobin_size=geobin_size, geobin_spacing=geobin_spacing,
+ margin=margin, filter=filter,
+ use_kde=True)
+ return d[1]*self.spacing[0], d[0]*self.spacing[1]
+
+ def shift_from_googlesat(self):
+ return
+
+ def apply_shift_inplace(self, shift):
+ self.scene = None
+ with rasterio.open(self.path, 'r+') as ds:
+ tr = ds.transform
+ with open(os.path.splitext(self.path)[0]+'_transform_backup.pkl', 'wb') as f:
+ pickle.dump(tr, f)
+ ds.transform = rasterio.transform.Affine(tr[0], tr[1], tr[2] + shift[0], tr[3], tr[4], tr[5] + shift[1])
+ self.scene = to_otb_pipeline(self.path)
+ self.spacing = self.scene.GetImageSpacing('out')
+ return