ENH: Sample points from S2 Planetary

TimeSeries/s2planetary.py

@@ -7,6 +7,10 @@ from Common.geotools import get_query_bbox
 from pyproj import Transformer as T
 from shapely.geometry import Polygon
 import tqdm
+from datetime import *
+import geopandas as gpd
+from scipy.interpolate import interp1d
+
 
 def fetch(shp, dt, output_fld, auth=None, band_list=None):
 
@@ -20,6 +24,7 @@ def fetch(shp, dt, output_fld, auth=None, band_list=None):
     if band_list is not None:
         lst = band_list
 
+    ofiles = []
     prg = tqdm.tqdm(total=len(res.item_collection()) * len(lst), desc="Fetching from Planetary")
     for item in res.items():
         with rasterio.open(item.assets['B02'].href) as ds:
@@ -34,6 +39,7 @@ def fetch(shp, dt, output_fld, auth=None, band_list=None):
 
         for a in lst:
             ofn = os.path.join(output_fld, '/'.join(item.assets[a].get_absolute_href().split('?')[0].split('/')[10:]))
+            ofiles.append(ofn)
             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,
@@ -54,7 +60,116 @@ def fetch(shp, dt, output_fld, auth=None, band_list=None):
         return S2PlaneteryPipeline(output_fld)
     else:
         warnings.warn("Queried for a non-default band list. Skipping pipeline setup.")
-        return
+        return ofiles
+
+def split_dt_monthly(dt):
+    sdt, edt = [datetime.strptime(x,'%Y-%m-%d') for x in dt.split('/')]
+    dts = [sdt]
+    cdt = sdt
+    while cdt < edt:
+        cdt += timedelta(weeks=4)
+        dts.append(cdt)
+    if dts[-1] > edt:
+        dts = dts[:-1]
+    if dts[-1] < edt:
+        dts.append(edt+timedelta(1))
+    out = []
+    for i in range(len(dts)-1):
+        out.append(datetime.strftime(dts[i], '%Y-%m-%d') + '/' + datetime.strftime(dts[i+1]-timedelta(1), '%Y-%m-%d'))
+    return out
+    
+def sample(shp, dt, band_list, field_name, temp_dir = '/tmp', strategy='all', auth = None):
+    fd = dt.split('/')[0]
+    fd = datetime.strptime(fd,'%Y-%m-%d')
+    ed = fd + timedelta(5)
+    ndt = datetime.strftime(fd, '%Y-%m-%d') + '/' + datetime.strftime(ed, '%Y-%m-%d')
+    #fn = fetch(shp, ndt, temp_dir, band_list=[band_list[0]])
+    fn = ['/work/temp/S2B_MSIL2A_20210104T103329_N0212_R108_T30PVT_20210121T134941.SAFE/GRANULE/L2A_T30PVT_A020011_20210104T104033/IMG_DATA/R10m/T30PVT_20210104T103329_B04_10m.tif',
+          '/work/temp/S2B_MSIL2A_20210104T103329_N0212_R108_T30PUT_20210105T133826.SAFE/GRANULE/L2A_T30PUT_A020011_20210104T104033/IMG_DATA/R10m/T30PUT_20210104T103329_B04_10m.tif']
+    
+    mos = otb.Registry.CreateApplication('Mosaic')
+    mos.SetParameterStringList('il', fn)
+    mos.Execute()
+    sts = otb.Registry.CreateApplication('PolygonClassStatistics')
+    sts.SetParameterInputImage('in', mos.GetParameterOutputImage('out'))
+    sts.SetParameterString('vec', shp)
+    sts.SetParameterString('field', field_name)
+    sts.SetParameterString('out', os.path.join(temp_dir, 'stats.xml'))
+    sts.ExecuteAndWriteOutput()
+
+    stats_fn = os.path.join(temp_dir, 'stats.xml')
+    samples_fn = os.path.join(temp_dir, 'samples.sqlite')
+
+    mos = otb.Registry.CreateApplication('Mosaic')
+    mos.SetParameterStringList('il', fn)
+    mos.Execute()
+    smp = otb.Registry.CreateApplication('SampleSelection')
+    smp.SetParameterInputImage('in', mos.GetParameterOutputImage('out'))
+    smp.SetParameterString('vec', shp)
+    smp.SetParameterString('instats', stats_fn)
+    smp.SetParameterString('field', field_name)
+    if isinstance(strategy, str):
+        smp.SetParameterString('strategy', strategy)
+    elif isinstance(strategy, int):
+        smp.SetParameterString('strategy','constant')
+        smp.SetParameterInt('strategy.constant.nb', strategy)
+    smp.SetParameterString('out', samples_fn)
+    smp.ExecuteAndWriteOutput()
+
+    ds = gpd.read_file(samples_fn)
+    coords = list(zip(list(ds.geometry.x),list(ds.geometry.y)))
+    bbox, _, _ = get_query_bbox(shp, return_all=True)
+    if auth is not None:
+        os.environ['PC_SDK_SUBSCRIPTION_KEY'] = auth
+    
+    all_arr = []
+    dates = []
+    for dti in split_dt_monthly(dt):
+
+        api = Client.open('https://planetarycomputer.microsoft.com/api/stac/v1', modifier=PC.sign_inplace)
+        res = api.search(collections="sentinel-2-l2a", bbox=bbox, datetime=dti, sortby='datetime')
+        found = res.item_collection()
+
+        if 'SCL' not in band_list:
+            band_list.append('SCL')
+        arr = {x : np.empty((ds.shape[0],len(found))) for x in band_list}
+        i = 0
+        prg = tqdm.tqdm(total=len(found), desc="Sampling from Planetary on period {}".format(dti))
+        for item in found:
+            dates.append(item.datetime)
+            for band in band_list:
+                with rasterio.open(item.assets[band].href) as src:
+                    arr[band][:,i] = np.array([x for x in src.sample(coords)]).flatten()
+            i += 1
+            prg.update()
+        prg.close()
+        api = None
+        all_arr.append(arr)
+
+    
+    days = np.array([(x-dates[0]).days for x in dates])
+    final_arr = {}
+    for b in band_list:
+        final_arr[b] = np.hstack([x[b] for x in all_arr])
+    
+    band_list.remove('SCL')
+    final_arr['SCL'] = np.invert(np.isin(final_arr['SCL'],[0,1,3,8,9,10]))
+    for i in tqdm.tqdm(range(ds.shape[0]), desc="Interpolating time series"):
+        x = days[final_arr['SCL'][i,:]]
+        for b in band_list:
+            y = final_arr[b][i,final_arr['SCL'][i,:]]
+            f = interp1d(x,y, fill_value='extrapolate')
+            final_arr[b][i,:] = f(days)
+    
+    _, ii = np.unique(days, return_index=True)
+    for b in band_list:
+        final_arr[b] = final_arr[b][:,ii]
+
+    #[os.remove(x) for x in fn]
+    os.remove(stats_fn)
+    os.remove(samples_fn)
+
+    return final_arr, dates[ii], list(ds[field_name.lower()])
 
 class S2PlanetaryTilePipeline(S2Sen2CorTilePipeline):
     NAME = 'S2-L2A-SEN2COR-PLANETARY'