Merge branch 'vector_bbox' into 'develop'

Vector bbox

See merge request remi.cresson/scenes!11

.gitignore

+.idea
+__pycache__

.gitlab-ci.yml

@@ -22,7 +22,7 @@ flake8:
 pylint:
   extends: .static_analysis_base
   script:
-    - sudo apt update && sudo apt install -y pylint && pylint --disable=too-many-nested-blocks,too-many-locals,too-many-statements,too-few-public-methods,too-many-instance-attributes,too-many-arguments,invalid-name,cell-var-from-loop,too-many-branches,too-many-ancestors --max-line-length=120 --ignored-modules=tqdm,rtree,pyotb $PWD/scenes
+    - sudo apt update && sudo apt install -y pylint && pylint --disable=too-many-locals,too-many-statements,too-few-public-methods,too-many-instance-attributes,invalid-name,too-many-branches --max-line-length=120 --ignored-modules=tqdm,rtree,pyotb $PWD/scenes
 
 codespell:
   extends: .static_analysis_base

scenes/drs.py

@@ -25,7 +25,7 @@ def get_spot67_scenes(root_dir):
     look_dir = root_dir + "/MS"
     print("Find files in {}".format(look_dir))
     dimap_xs_files = find_all_dimaps(look_dir)
-    print("Found {} DIMAP files in MS folder".format(len(dimap_xs_files)))
+    print("Found {} DIMAP files in MS folder".format(len(dimap_xs_files)), flush=True)
 
     # Create scenes list
     scenes = []

scenes/utils.py

@@ -4,6 +4,7 @@ A set of helpers for generic purpose
 import os
 import glob
 import pathlib
+import math
 from osgeo import osr, ogr, gdal
 
 
@@ -80,17 +81,25 @@ def extent_to_bbox(extent):
     :param extent: extent
     :return: bounding box (xmin, xmax, ymin, ymax)
     """
-    (xmin, ymax), (xmax, _), (_, ymin), (_, _) = extent
-    return min(xmin, xmax), max(xmin, xmax), min(ymin, ymax), max(ymin, ymax)
+    xmin, xmax = math.inf, -math.inf
+    ymin, ymax = math.inf, -math.inf
+    for x, y in extent:
+        xmin = min(xmin, x)
+        xmax = max(xmax, x)
+        ymin = min(ymin, y)
+        ymax = max(ymax, y)
 
+    return xmin, xmax, ymin, ymax
 
-def get_bbox_wgs84(gdal_ds):
+
+def get_bbox_wgs84_from_gdal_ds(gdal_ds):
     """
     Returns the bounding box in WGS84 CRS from a GDAL dataset
     :param gdal_ds: GDAL dataset
     :return: bounding box (xmin, xmax, ymin, ymax)
     """
     extend_wgs84 = get_extent_wgs84(gdal_ds)
+
     return extent_to_bbox(extend_wgs84)
 
 
@@ -104,6 +113,7 @@ def get_epsg_extent_bbox(filename):
     epsg = get_epsg(gdal_ds)
     extent_wgs84 = get_extent_wgs84(gdal_ds)
     bbox_wgs84 = extent_to_bbox(extent_wgs84)
+
     return epsg, extent_wgs84, bbox_wgs84
 
 
@@ -166,16 +176,32 @@ def extent_overlap(extent, other_extent):
 
 def open_vector_layer(vector_file):
     """
-    Return the vector data layer from a vector file. If the vector is empty, None is returned.
+    Return the vector dataset from a vector file. If the vector is empty, None is returned.
     :param vector_file: input vector file
-    :return: ogr layer, or None (if layer is empty)
+    :return: ogr ds, or None (if error)
     """
     poly_ds = ogr.Open(vector_file)
     if poly_ds is None:
         raise Exception("ERROR: vector layer {} is None!".format(vector_file))
-    poly_layer = poly_ds.GetLayer()
+    return poly_ds
+
 
-    return poly_layer
+def get_bbox_wgs84_from_vector(vector_file):
+    """
+    Returns the bounding box in WGS84 CRS from a vector data
+    :param vector_file: vector data filename
+    :return: bounding box in WGS84 CRS
+    """
+    poly_ds = open_vector_layer(vector_file=vector_file)
+    poly_layer = poly_ds.GetLayer()
+    extent = poly_layer.GetExtent()
+    coords = [(extent[0], extent[2]), (extent[1], extent[3])]
+    src_srs = poly_layer.GetSpatialRef()
+    tgt_srs = epsg2srs(4326)
+    [(xmin, ymin), (xmax, ymax)] = reproject_coords(coords=coords,  # pylint: disable=unbalanced-tuple-unpacking
+                                                    src_srs=src_srs,
+                                                    tgt_srs=tgt_srs)
+    return xmin, xmax, ymin, ymax
 
 
 def find_file_in_all_subdirs(pth, pattern):

test/drs_search.py

@@ -18,8 +18,11 @@ idx = indexation.Index(scenes)
 
 # search
 print("search roi")
-gdal_ds = gdal.Open(params.roi)
-bbox = utils.get_bbox_wgs84(gdal_ds=gdal_ds)
+if params.roi.lower().endswith(".tif"):
+    gdal_ds = gdal.Open(params.roi)
+    bbox = utils.get_bbox_wgs84_from_gdal_ds(gdal_ds=gdal_ds)
+else:
+    bbox = utils.get_bbox_wgs84_from_vector(params.roi)
 matches = idx.find(bbox_wgs84=bbox)
 print("{} scenes found.".format(len(matches)))
 #for scene_match in matches:

test/indexation_test.py

 # -*- coding: utf-8 -*-
 from scenes_test_base import ScenesTestBase
-from scenes import indexation
+from scenes import indexation, utils
 import tests_data
 
 class ImageryTest(ScenesTestBase):
@@ -9,6 +9,8 @@ class ImageryTest(ScenesTestBase):
         index = indexation.Index(scenes_list=[tests_data.SCENE1])
         self.assertTrue(index.find(bbox_wgs84=(43.706, 43.708, 4.317, 4.420)))
         self.assertFalse(index.find(bbox_wgs84=(43.000, 43.001, 3.000, 3.001)))
+        self.assertTrue(index.find(bbox_wgs84=utils.get_bbox_wgs84_from_vector(tests_data.ROI_MTP_4326)))
+        self.assertTrue(index.find(bbox_wgs84=utils.get_bbox_wgs84_from_vector(tests_data.ROI_MTP_2154)))
 
     def test_epsg(self):
         self.assertTrue(tests_data.SCENE1.epsg == 2154)

test/tests_data.py

@@ -11,6 +11,8 @@ DIMAP1_XS = TEST_DATA_DIR + "/input/ROI_1_Bundle_Ortho_GSD2015/PROD_SPOT6_001/VO
                             "DIM_SPOT6_MS_201503261014386_ORT_SPOT6_20170524_1422391k0ha487979cy_1.XML"
 DIMAP1_P = TEST_DATA_DIR + "/input/ROI_1_Bundle_Ortho_GSD2015/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_P_001_A/" \
                            "DIM_SPOT6_P_201503261014386_ORT_SPOT6_20170524_1422391k0ha487979cy_1.XML"
+ROI_MTP_2154 = TEST_DATA_DIR + "/input/roi_mtp_2154.gpkg"
+ROI_MTP_4326 = TEST_DATA_DIR + "/input/roi_mtp_4326.gpkg"
 
 # Instances
 SCENE1 = spot.Spot67Scene(dimap_file_xs=DIMAP1_XS, dimap_file_pan=DIMAP1_P)
\ No newline at end of file