diff --git a/TimeSeries/s1base.py b/TimeSeries/s1base.py
index 3a53d80936f37c59660546bb6e045536a4abf95e..d0cbd0977f20ecc56b809f40338e67e3eecc813f 100644
--- a/TimeSeries/s1base.py
+++ b/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:
diff --git a/TimeSeries/s1planetary.py b/TimeSeries/s1planetary.py
new file mode 100644
index 0000000000000000000000000000000000000000..e104c9303c0fcb4567c183525a81c6e229b21ae8
--- /dev/null
+++ b/TimeSeries/s1planetary.py
@@ -0,0 +1,108 @@
+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