ENH: Can mine folders containing tiles non overlapping with ROI.

Common/geotools.py

 from osgeo import gdal
 from osgeo import ogr
 from pyproj import Transformer as T
+from shapely.geometry import box
+from shapely.ops import transform
+import rasterio
+import geopandas
 
 def get_query_bbox(shp, query_srs=4326, return_all=False):
     ds = ogr.Open(shp)
@@ -20,4 +24,19 @@ def get_query_bbox(shp, query_srs=4326, return_all=False):
     if return_all:
         return bbox, i_bbox, shp_srs
     else:
-        return bbox
\ No newline at end of file
+        return bbox
+    
+def check_extent_overlap(im,vect):
+    im_ds = rasterio.open(im)
+    im_srs = im_ds.crs
+    im_extent = box(*im_ds.bounds)
+
+    vect_ds = geopandas.read_file(vect)
+    vect_srs = vect_ds.crs
+    vect_extent = box(*vect_ds.geometry.total_bounds)
+
+    if (im_srs != vect_srs):
+        tr = T.from_crs(im_srs, vect_srs, always_xy=True)
+        im_extent = transform(tr.transform, im_extent)
+        
+    return im_extent.intersects(vect_extent)

TimeSeries/s2theia.py

@@ -14,7 +14,7 @@ from osgeo import osr
 import datetime
 import uuid
 import shutil
-from Common.geotools import get_query_bbox
+from Common.geotools import get_query_bbox, check_extent_overlap
 from Common.geometry import get_displacements_to_ref, get_clearest_central_image
 
 def fetch(shp, dt, output_fld, credentials):
@@ -145,9 +145,8 @@ class S2TheiaTilePipeline:
 
     @classmethod
     def _check_roi(cls, x, roi, min_surf=0.0, temp_fld='/tmp'):
-        if cls._check(x):
+        if cls._check(x):            
             if os.path.isdir(x):
-                er = otb.Registry.CreateApplication('ExtractROI')
                 bnd = glob.glob(os.path.join(x, cls.PTRN_20m[0]))[0]
             elif os.path.splitext(x)[-1] == '.zip':
                 idf = cls.PTRN_20m[0].replace('*', '')
@@ -158,6 +157,9 @@ class S2TheiaTilePipeline:
                 tgt = open(bnd, 'wb')
                 with fle,tgt:
                     shutil.copyfileobj(fle, tgt)
+            if not check_extent_overlap(bnd, roi):
+                return False
+            er = otb.Registry.CreateApplication('ExtractROI')
             er.SetParameterString('in', bnd)
             er.SetParameterString('mode', 'fit')
             er.SetParameterString('mode.fit.vect', roi)
@@ -219,24 +221,30 @@ class S2TheiaTilePipeline:
                     self.input_dates = [self.input_dates[i] for i in idx]
                     self.tile_cloud_percentage = [self.tile_cloud_percentage[i] for i in idx]
 
-                self.set_input_epsg()
-                self.output_epsg = self.input_epsg
-                self.output_dates = self.input_dates
-
-                for img in self.image_list:
-                    for p in self.PTRN_ful:
-                        ifn, ofn = self.get_file(img, p)
-                        self.append(to_otb_pipeline(ifn), ofn, self.REF_TYPE, 'out', is_output=True)
-                    for p in self.PTRN_msk:
-                        ifn, ofn = self.get_file(img, p)
-                        self.append(to_otb_pipeline(ifn), ofn, self.MSK_TYPE, 'out', is_output=True)
+                if len(self.image_list) > 0:
+                    self.set_input_epsg()
+                    self.output_epsg = self.input_epsg
+                    self.output_dates = self.input_dates
+
+                    for img in self.image_list:
+                        for p in self.PTRN_ful:
+                            ifn, ofn = self.get_file(img, p)
+                            self.append(to_otb_pipeline(ifn), ofn, self.REF_TYPE, 'out', is_output=True)
+                        for p in self.PTRN_msk:
+                            ifn, ofn = self.get_file(img, p)
+                            self.append(to_otb_pipeline(ifn), ofn, self.MSK_TYPE, 'out', is_output=True)
+                else:
+                    warnings.warn('Empty pipeline after filtering for tile {}.'.format(self.tile_id))
 
         else:
-            warnings.warn('Empty pipeline. Need to set preprocessed inputs?')
+            warnings.warn('Empty pipeline at init. Need to set preprocessed inputs?')
 
     def __del__(self):
         if os.path.exists(self.temp_fld):
             shutil.rmtree(self.temp_fld)
+    
+    def is_empty(self):
+        return len(self.image_list) == 0
 
     def reset(self):
         self.pipe = []
@@ -775,20 +783,25 @@ class S2TheiaPipeline:
     _tile_id = S2TilePipeline._tile_id
     tiles = []
 
-    def __init__(self, fld, temp_fld='/tmp', input_date_interval=None, max_clouds_percentage=None):
+    def __init__(self, fld, temp_fld='/tmp', input_date_interval=None, max_clouds_percentage=None, roi=None):
         self.folder = fld
         self.temp_fld = temp_fld
         self.input_date_interval = input_date_interval
         self.max_clouds_percentage = max_clouds_percentage
         self.tile_list = set()
+        self.roi = roi
         img_list = [os.path.abspath(x) for x in glob.glob(os.path.join(self.folder, self.S2TilePipeline.PTRN_dir))
                     if os.path.isdir(x) or os.path.splitext(x)[-1] == '.zip']
         [self.tile_list.add(self._tile_id(x)) for x in img_list if self._check(x)]
         assert all([self.S2TilePipeline(fld, t, dummy_read=True) for t in self.tile_list])
-        self.tiles = [self.S2TilePipeline(fld, t, self.temp_fld, self.input_date_interval, self.max_clouds_percentage) for t in self.tile_list]
+        self.tiles = [self.S2TilePipeline(fld, t, self.temp_fld, self.input_date_interval, self.max_clouds_percentage, filter_by_roi=self.roi) for t in self.tile_list]
+        self.tiles = [x for x in self.tiles if not x.is_empty()]
         self.output_dates = None
-        self.roi = None
-        self.output_epsg = self.tiles[0].input_epsg
+
+        if len(self.tiles) > 0:
+            self.output_epsg = self.tiles[0].input_epsg
+        else:
+            warnings.warn('Empty pipeline.')
 
     def __del__(self):
         for x in self.tiles:

Workflows/operations.py

@@ -40,9 +40,9 @@ def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None,
     else:
         raise ValueError("Unsupported/non-valid provider")
 
-    s2 = S2Processor(in_fld, temp_fld=out_fld)
-    if roi is not None:
-        s2.set_roi(roi)
+    s2 = S2Processor(in_fld, temp_fld=out_fld, roi=roi)
+    #if roi is not None:
+    #    s2.set_roi(roi)
     if output_dates_file is not None:
         s2.set_output_dates_by_file(output_dates_file)
     else: