diff --git a/Common/geotools.py b/Common/geotools.py
new file mode 100644
index 0000000000000000000000000000000000000000..646c2d2db27795f091f71e784baf0a9defcf8526
--- /dev/null
+++ b/Common/geotools.py
@@ -0,0 +1,23 @@
+from osgeo import gdal
+from osgeo import ogr
+from pyproj import Transformer as T
+
+def get_query_bbox(shp, query_srs=4326, return_all=False):
+ ds = ogr.Open(shp)
+ ly = ds.GetLayer(0)
+ xt = ly.GetExtent()
+ shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))
+
+ i_bbox = [xt[0], xt[2], xt[1], xt[3]]
+ if shp_srs == query_srs:
+ bbox = i_bbox
+ else:
+ tr = T.from_crs(shp_srs, query_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, lat_min, lon_max, lat_max]
+
+ if return_all:
+ return bbox, i_bbox, shp_srs
+ else:
+ return bbox
\ No newline at end of file
diff --git a/TimeSeries/s1base.py b/TimeSeries/s1base.py
index d0cbd0977f20ecc56b809f40338e67e3eecc813f..923b2f0d8eeb741340599eb3ba5d74ed2cd3b0af 100644
--- a/TimeSeries/s1base.py
+++ b/TimeSeries/s1base.py
@@ -7,13 +7,34 @@ import uuid
import xml.etree.ElementTree as ET
import shutil
from itertools import groupby
+from Common.geotools import get_query_bbox
+from eodag import EODataAccessGateway
+
+def fetch(shp, dt, output_fld, credentials):
+ bbox = get_query_bbox(shp)
+ dag = EODataAccessGateway(user_conf_file_path=credentials)
+ dag.set_preferred_provider("peps")
+ search_criteria = {
+ "productType": "S1_SAR_GRD",
+ "start": dt.split("/")[0],
+ "end": dt.split("/")[1],
+ "geom": {"lonmin": bbox[0], "latmin": bbox[1], "lonmax": bbox[2], "latmax": bbox[3]}
+ }
+ res = dag.search_all(**search_criteria)
+ res.filter_property(sensorMode='IW')
+ if len(res) > 0:
+ os.makedirs(output_fld, exist_ok=True)
+ dag.download_all(res, outputs_prefix=output_fld, extract=True)
+
+ return S1GRDPipeline(output_fld)
+
class S1GRDPipeline:
# --- BEGIN SENSOR PROTOTYPE ---
NAME = 'S1-IW-GRD'
VAL_TYPE = otb.ImagePixelType_int16
TMP_TYPE = otb.ImagePixelType_float
- PTRN_dir = 'S1*_IW_GRD*.SAFE'
+ PTRN_dir = 'S1*_IW_GRD*/S1*_IW_GRD*.SAFE'
PTRN_ref = '-iw-grd-'
VH_name = 'vh'
PTRN = ['measurement/s1*-iw-grd-vh-*.tiff', 'measurement/s1*-iw-grd-vv-*.tiff']
diff --git a/TimeSeries/s1planetary.py b/TimeSeries/s1planetary.py
index e104c9303c0fcb4567c183525a81c6e229b21ae8..28e04c2685f4fd36684e5ef2a3a85b013e3eb0eb 100644
--- a/TimeSeries/s1planetary.py
+++ b/TimeSeries/s1planetary.py
@@ -1,9 +1,7 @@
import os.path
-
from TimeSeries.s1base import *
import planetary_computer as PC
from pystac_client import Client
-from osgeo import ogr
from pyproj import Transformer as T
from shapely.geometry import Polygon
import rasterio
@@ -11,24 +9,12 @@ import rasterio.mask
import tqdm
import otbApplication as otb
import json
+from Common.geotools import get_query_bbox
def fetch(shp, dt, output_fld, api_key, direction=None):
os.environ['PC_SDK_SUBSCRIPTION_KEY'] = api_key
- ds = ogr.Open(shp)
- ly = ds.GetLayer(0)
- xt = ly.GetExtent()
- shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))
- qry_srs = 4326
-
- i_bbox = [xt[0], xt[2], xt[1], xt[3]]
- if shp_srs == qry_srs:
- bbox = i_bbox
- else:
- tr = T.from_crs(shp_srs, qry_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, lat_min, lon_max, lat_max]
+ bbox, i_bbox, shp_srs = get_query_bbox(shp, return_all=True)
api = Client.open('https://planetarycomputer.microsoft.com/api/stac/v1', modifier=PC.sign_inplace)
res = api.search(collections="sentinel-1-rtc", bbox=bbox, datetime=dt)
diff --git a/TimeSeries/s2planetary.py b/TimeSeries/s2planetary.py
index a0de40366d7655da645181a4762821b2376245fa..928ce5399ea04c0120d91aec304639f6924d6e33 100644
--- a/TimeSeries/s2planetary.py
+++ b/TimeSeries/s2planetary.py
@@ -1,31 +1,16 @@
-import warnings
-
from TimeSeries.s2sen2cor import *
import planetary_computer as PC
from pystac_client import Client
import rasterio
import rasterio.mask
-from osgeo import ogr
+from Common.geotools import get_query_bbox
from pyproj import Transformer as T
from shapely.geometry import Polygon
import tqdm
def fetch(shp, dt, output_fld, band_list=None):
- ds = ogr.Open(shp)
- ly = ds.GetLayer(0)
- xt = ly.GetExtent()
- shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))
- qry_srs = 4326
-
- i_bbox = [xt[0], xt[2], xt[1], xt[3]]
- if shp_srs == qry_srs:
- bbox = i_bbox
- else:
- tr = T.from_crs(shp_srs,qry_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, lat_min, lon_max, lat_max]
+ bbox, i_bbox, shp_srs = get_query_bbox(shp, return_all=True)
api = Client.open('https://planetarycomputer.microsoft.com/api/stac/v1', modifier=PC.sign_inplace)
res = api.search(collections="sentinel-2-l2a", bbox=bbox, datetime=dt)
diff --git a/TimeSeries/s2theia.py b/TimeSeries/s2theia.py
index 39593e2c02b22f55eae467783853a3f5bbdeecc6..bd9ff9f0b2035064dbc3eb296f1fb2426d5ecfce 100644
--- a/TimeSeries/s2theia.py
+++ b/TimeSeries/s2theia.py
@@ -1,11 +1,7 @@
-import sys
import warnings
-
-from osgeo import gdal,ogr
+from osgeo import gdal
import otbApplication as otb
from theia_picker import TheiaCatalog
-from pyproj import Transformer as T
-import tqdm
from Common.otb_numpy_proc import to_otb_pipeline
import numpy as np
@@ -17,25 +13,11 @@ from osgeo import osr
import datetime
import uuid
import shutil
-
+from Common.geotools import get_query_bbox
from Common.geometry import get_displacements_to_ref, get_clearest_central_image
def fetch(shp, dt, output_fld, credentials):
-
- ds = ogr.Open(shp)
- ly = ds.GetLayer(0)
- xt = ly.GetExtent()
- shp_srs = int(ly.GetSpatialRef().GetAuthorityCode(None))
- qry_srs = 4326
-
- i_bbox = [xt[0], xt[2], xt[1], xt[3]]
- if shp_srs == qry_srs:
- bbox = i_bbox
- else:
- tr = T.from_crs(shp_srs,qry_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, lat_min, lon_max, lat_max]
+ bbox = get_query_bbox(shp)
theia = TheiaCatalog(credentials)
features = theia.search(
@@ -48,11 +30,8 @@ def fetch(shp, dt, output_fld, credentials):
lst = ['FRE_B2','FRE_B3','FRE_B4','FRE_B5','FRE_B6','FRE_B7','FRE_B8','FRE_B8A',
'FRE_B11','FRE_B12','EDG_R1','SAT_R1','CLM_R1']
- #prg = tqdm.tqdm(total=len(features) * 13, desc="Fetching from Theia")
for f in features:
f.download_files(matching=lst, download_dir=output_fld)
- #prg.update()
- #prg.close()
return S2TheiaPipeline(output_fld)
diff --git a/docker/dockerfile b/docker/dockerfile
index 38766075c4abfa6bab8bdef5beba6a8453ce4cad..47c533127a544f19c49781e5ee1ced33e1e492a4 100644
--- a/docker/dockerfile
+++ b/docker/dockerfile
@@ -115,7 +115,8 @@ RUN pip install gdal==3.4.2 \
planetary_computer \
theia_picker \
fpdf2 \
- planet
+ planet \
+ eodag
RUN mkdir moringav2
COPY --chown=ubuntu . /home/ubuntu/moringav2
\ No newline at end of file
diff --git a/moringa.py b/moringa.py
index f297185e99d2272e46458a4ebd8cde05ad7f1066..c372a65449a98a8c60ee7630d4b5c64d593636f7 100755
--- a/moringa.py
+++ b/moringa.py
@@ -3,7 +3,7 @@ import sys
import argparse
import OBIA.segmentation
import VHR.vhrbase
-from TimeSeries import s2theia, s2planetary
+from TimeSeries import s2theia, s2planetary, s1base, s1planetary
def run_segmentation(img, threshold, cw, sw , out_seg,
n_first_iter, margin, roi, n_proc, memory,
@@ -32,8 +32,8 @@ def preprocess_spot67(in_fld, out_fld, dem_fld, geoid_file, skip_ps, compress,
sp.write_outputs(out_fld, compress=compress)
return
-def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None, coregister_to=None, coregister_to_band=1,
- coregister_using_band=3, provider='theia'):
+def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None,
+ align_to=None, align_to_band=3, align_using_band=3, provider='theia'):
S2Processor = None
if provider == 'theia':
S2Processor = s2theia.S2TheiaPipeline
@@ -50,23 +50,50 @@ def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None, coregister_
else:
raise ValueError("Please provide path to a text file containing output dates.")
- s2.extract_feature_set(out_fld, store_gapfill=True, mosaicking='vrt')
+ align = align_to is not None
+ if (align_to == 'self'):
+ align_to = None
+
+ s2.extract_feature_set(out_fld, store_gapfill=True, mosaicking='vrt', align=align,
+ align_to=align_to, align_to_band=align_to_band, align_using_band=align_using_band)
+ return
+
+def preprocess_s1(in_fld, roi, out_fld, dem_fld=None, geoid=None, direction=None, satellite=None,
+ skip_despeckle=False, provider='native'):
+ S1processor = None
+ if provider == 'native':
+ S1Processor = s1base.S1GRDPipeline
+ elif provider == 'planetary':
+ S1Processor = s1planetary.S1RTCPipeline
+ else:
+ raise ValueError("Unsupported/non-valid provider")
+
+ s1 = S1Processor(in_fld, roi, out_fld, direction=direction, satellite=satellite)
+ if provider == 'native':
+ assert(os.path.exists(dem_fld) and os.path.exists(geoid))
+ s1.calibrate()
+ s1.orthorectify(dem_fld, geoid)
+ s1.stitch()
+ if not skip_despeckle:
+ s1.multitemp_speckle_filter()
+ s1.compute_features()
+ s1.write_outputs(out_fld)
return
def main(args):
parser = argparse.ArgumentParser(prog="moringa", add_help=False)
subpar = parser.add_subparsers(dest="cmd")
- prepr = subpar.add_parser("preprocess_s2", help="Performs Moringa preset time series preprocessing.",
+ prepr = subpar.add_parser("preprocess_s2", help="Performs Moringa preset time series preprocessing for Sentinel-2.",
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
prepr.add_argument("in_folder", type=str, help="Path to the folder containing (de-zipped) S2 images.")
prepr.add_argument("out_folder", type=str, help="Path to the folder in which preprocessed stacks will be stored.")
prepr.add_argument("--output_dates_file", type=str, default=None, help="Path to the text file containing output dates for temporal interpolation.")
prepr.add_argument("--roi", type=str, default=None, help="Path to the ROI vector file.")
- prepr.add_argument("--coregister_to", type=str, default=None, help="Path to a reference image to which the stacks must be coregistered.")
- prepr.add_argument("--coregister_to_band", type=int, default=1, help="Band of reference image used for co-registration.")
- prepr.add_argument("--coregister_using_band", type=int, default=3, help="Band of current stack used for co-registration.")
- prepr.add_argument("--provider", type=str, default='theia', help="S2 image provider. Supported: 'theia', 'theial3a', 'sen2cor', 'planetary'")
+ prepr.add_argument("--align_to", type=str, default=None, help="Optional strategy for S2 spatial alignment (self/<date>/path to reference image).")
+ prepr.add_argument("--align_to_band", type=int, default=3, help="Band of reference image used for co-registration.")
+ prepr.add_argument("--align_using_band", type=int, default=3, help="Band of current stack used for co-registration.")
+ prepr.add_argument("--provider", type=str, default='theia', help="S2 image provider. Currently supported: 'theia', 'planetary'")
segmt = subpar.add_parser("segment", help="Performs (large scale Baatz-Shape) segmentation of an input image.",
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
@@ -97,7 +124,18 @@ def main(args):
vhrprep.add_argument("--skip_ps", help="Skip pansharpening step.", action='store_true')
vhrprep.add_argument("--compress", help="Use lossless compression on outputs.", action='store_true')
-
+ s1prepr = subpar.add_parser("preprocess_s1", help="Performs Moringa preset time series preprocessing for Sentinel-1.",
+ formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+ s1prepr.add_argument("in_folder", type=str, help="Path to the folder containing (de-zipped) S1 images.")
+ s1prepr.add_argument("roi", type=str, default=None, help="Path to an image whose geometry has to be fit (e.g. a S2 scene).")
+ s1prepr.add_argument("out_folder", type=str, help="Path to the folder in which preprocessed stacks will be stored.")
+ s1prepr.add_argument("--dem_fld", type=str, default=None, help="Path to the folder containing DEM covering the scene in WGS84 projection.")
+ s1prepr.add_argument("--geoid", type=str, default=None, help="Path to the geoid file.")
+ s1prepr.add_argument("--direction", type=str, default=None, help="Filter direction (ascending/descending)")
+ s1prepr.add_argument("--satellite", type=str, default=None, help="Filter satellite (s1a/s1b)")
+ s1prepr.add_argument("--skip_despeckle", help="Skip despeckling step", action='store_true')
+ s1prepr.add_argument("--provider", type=str, default='native',
+ help="S1 image provider. Currently supported: 'native' (e.g. esa/peps), 'planetary'")
if len(args) == 1:
parser.print_help()
@@ -107,7 +145,8 @@ def main(args):
if arg.cmd == "preprocess_s2":
preprocess_s2(arg.in_folder, arg.out_folder, output_dates_file=arg.output_dates_file, roi=arg.roi,
- coregister_to=arg.coregister_to, provider=arg.provider)
+ align_to=arg.align_to, align_to_band=arg.align_to_band, align_using_band=arg.align_using_band,
+ provider=arg.provider)
if arg.cmd == "segment":
run_segmentation(arg.img, arg.threshold, arg.cw, arg.sw, arg.outimg, arg.n_first_iter, arg.tile_margin,
@@ -117,6 +156,10 @@ def main(args):
preprocess_spot67(arg.fld, arg.out_fld, arg.dem_fld, arg.geoid, arg.skip_ps, arg.compress,
arg.clip, arg.align_to, arg.align_to_band, arg.align_using_band)
+ if arg.cmd == "preprocess_s1":
+ preprocess_s1(arg.fld, arg.roi, arg.out_fld, arg.dem_fld, arg.geoid, arg.direction, arg.satellite,
+ arg.skip_despeckle, arg.provider)
+
return 0
if __name__ == "__main__":