ENH: improve S1 pipeline and integrate Planetary RTC product.

TimeSeries/s1base.py

@@ -34,14 +34,14 @@ class S1GRDPipeline:
     @classmethod
     def _img_id(cls, x):
         if cls._check(x):
-            return '_'.join(x.split('_')[4:7])
+            return '_'.join(os.path.basename(x).split('_')[4:7])
         else:
             return None
 
     @classmethod
     def _img_date(cls, x):
         if cls._check(x):
-            return x.split('_')[4].split('T')[0]
+            return os.path.basename(x).split('_')[4].split('T')[0]
         else:
             return None
 
@@ -76,10 +76,16 @@ class S1GRDPipeline:
                 else:
                     return None
 
+    @classmethod
+    def _get_stiched_filename(cls, fn):
+        tm1 = os.path.basename(fn).split('-')[-4].split('t')[1]
+        tm2 = os.path.basename(fn).split('-')[-5].split('t')[1]
+        fn = fn.replace(tm1, 'xxxxxx').replace(tm2, 'xxxxxx')
+        return fn
+
     def _check_satellite(self, x):
         return os.path.basename(x).split('_')[0].lower()
 
-
     def __init__(self, fld, roi, temp_fld='/tmp', input_date_interval=None, roi_min_surf=0.0,
                  direction=None, satellite=None):
         self.pipe = []
@@ -198,10 +204,7 @@ class S1GRDPipeline:
                             self.types[proc_idx[2 * s[0] + 1]], self.out_p[proc_idx[2 * s[0] + 1]], is_output=True)
             else:
                 for k in range(2):
-                    fn = self.files[proc_idx[2*s[0]+k]]
-                    tm1 = os.path.basename(fn).split('-')[-4].split('t')[1]
-                    tm2 = os.path.basename(fn).split('-')[-5].split('t')[1]
-                    fn = fn.replace(tm1, 'xxxxxx').replace(tm2, 'xxxxxx')
+                    fn = self._get_stitched_filename(self.files[proc_idx[2*s[0]+k]])
                     bm = otb.Registry.CreateApplication('BandMathX')
                     [bm.AddImageToParameterInputImageList('il', self.pipe[proc_idx[2*u+k]].GetParameterOutputImage(self.out_p[proc_idx[2*u+k]])) for u in s]
                     bm.SetParameterString('exp',
@@ -268,7 +271,7 @@ class S1GRDPipeline:
             expr = '{' + ';'.join(expr) + '}'
             bm.SetParameterString('exp', expr)
             bm.Execute()
-            fn = self.files[proc_idx[k]].replace('-vh-', '-feat-')
+            fn = self.files[proc_idx[k]].replace('-vh', '-feat')
             if all(['db' in x for x in feat_list]):
                 ty = self.VAL_TYPE
             else:

TimeSeries/s1planetary.py

+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
+import rasterio.mask
+import tqdm
+import otbApplication as otb
+import json
+
+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]
+
+    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)
+    lst = ['vh', 'vv']
+
+    prg = tqdm.tqdm(total=len(res.item_collection()) * len(lst), desc="Fetching from Planetary")
+    for item in res.items():
+        if direction is not None:
+            if item.properties['sat:orbit_state'] != direction:
+                continue
+        mfn = os.path.join(output_fld, item.assets['vh'].get_absolute_href().split('?')[0].split('/')[10])
+        mfn = os.path.join(mfn, 'item_properties.json')
+        os.makedirs(os.path.dirname(mfn), exist_ok=True)
+        with open(mfn, 'w') as js:
+            json.dump(item.properties, js)
+        with rasterio.open(item.assets['vh'].href) as ds:
+            img_srs = ds.crs.to_epsg()
+        if shp_srs == img_srs:
+            o_bbox = [20 * round(x / 20) for x in i_bbox]
+        else:
+            tr = T.from_crs(shp_srs, img_srs, always_xy=True)
+            x_min, y_min = tr.transform(i_bbox[0], i_bbox[1])
+            x_max, y_max = tr.transform(i_bbox[2], i_bbox[3])
+            o_bbox = [20 * round(x / 20) for x in [x_min, y_min, x_max, y_max]]
+
+        for a in lst:
+            ofn = os.path.join(output_fld, '/'.join(item.assets[a].get_absolute_href().split('?')[0].split('/')[10:]))
+            os.makedirs(os.path.dirname(ofn), exist_ok=True)
+            with rasterio.open(item.assets[a].href) as img:
+                out_img, out_geot = rasterio.mask.mask(img,
+                                                       [Polygon.from_bounds(o_bbox[0], o_bbox[1],
+                                                                            o_bbox[2], o_bbox[3])],
+                                                       crop=True)
+                out_meta = img.meta
+            out_meta.update({"driver": "GTiff",
+                             "height": out_img.shape[1],
+                             "width": out_img.shape[2],
+                             "transform": out_geot})
+            with rasterio.open(ofn, "w", **out_meta) as dest:
+                dest.write(out_img)
+            prg.update()
+    prg.close()
+
+    return S1RTCPipeline(output_fld, ofn, direction=direction)
+
+class S1RTCPipeline(S1GRDPipeline):
+    # --- BEGIN SENSOR PROTOTYPE ---
+    NAME = 'S1-IW-RTC'
+    VAL_TYPE = otb.ImagePixelType_int16
+    PTRN_dir = 'S1*_IW_GRD*'
+    PTRN_ref = 'iw-v'
+    VH_name = 'vh'
+    PTRN = ['measurement/iw-vh.rtc.tiff', 'measurement/iw-vv.rtc.tiff']
+    FEAT_exp = {
+        'VH': 'im1b1',
+        'VV': 'im2b1',
+        'VH_db': '1000*log10(abs(im1b1)+1e-6)',
+        'VV_db': '1000*log10(abs(im2b1)+1e-6)',
+        'POL_RATIO': 'im1b1/im2b1'
+    }
+    NDT = -32768.0
+
+    @classmethod
+    def _check_direction(cls, x):
+        with open(os.path.join(x, 'item_properties.json'), 'r') as js:
+            dct = json.load(js)
+        if dct['sat:orbit_state'] == 'ascending':
+            return 0
+        elif dct['sat:orbit_state'] == 'descending':
+            return 1
+        else:
+            return None
+
+    @classmethod
+    def _get_stitched_filename(cls, fn):
+        fnx = os.path.splitext(fn)
+        return fnx[0] + '.stitched' + fnx[1]
\ No newline at end of file