Commit 1eb14525 authored by Cresson Remi's avatar Cresson Remi
Browse files

Merge branch '13-bbox_from_otbobj' into 'develop'

BoundingBox from otbObjects, GDAL dataset, or filename

Closes #13

See merge request !31
1 merge request!31BoundingBox from otbObjects, GDAL dataset, or filename
Pipeline #35094 passed with stages
in 4 minutes and 22 seconds
Showing with 76 additions and 15 deletions
+76 -15
...@@ -6,7 +6,6 @@ from scenes.spot import get_spot67_scenes ...@@ -6,7 +6,6 @@ from scenes.spot import get_spot67_scenes
# Arguments # Arguments
parser = argparse.ArgumentParser(description="Import Spot-6/7 images from DRS into a list of scenes (pickle)",) parser = argparse.ArgumentParser(description="Import Spot-6/7 images from DRS into a list of scenes (pickle)",)
parser.add_argument("--root_dirs", nargs='+', help="Root directories containing MS and PAN folders", required=True) parser.add_argument("--root_dirs", nargs='+', help="Root directories containing MS and PAN folders", required=True)
parser.add_argument("--tile_name", help="Filter only images from this tile (optional)")
parser.add_argument("--out_pickle", help="Output pickle file", required=True) parser.add_argument("--out_pickle", help="Output pickle file", required=True)
params = parser.parse_args() params = parser.parse_args()
......
""" """
This module contains a set of functions to deal with GDAL datasets (rasters) This module contains a set of functions to deal with GDAL datasets (rasters)
""" """
import pyotb
from osgeo import osr, gdal from osgeo import osr, gdal
from scenes.vector import reproject_coords from scenes.vector import reproject_coords
from scenes.spatial import epsg2srs, coord_list_to_bbox from scenes.spatial import epsg2srs, coord_list_to_bbox
...@@ -23,39 +24,69 @@ def get_epsg(gdal_ds): ...@@ -23,39 +24,69 @@ def get_epsg(gdal_ds):
return int(epsg) return int(epsg)
def get_extent(gdal_ds): def get_extent(raster):
""" """
Return list of corners coordinates from a GDAL dataset Return list of corners coordinates from a raster
Args: Args:
gdal_ds: GDAL dataset raster: GDAL dataset, str (raster filename) or otbobject
Returns: Returns:
list of coordinates list of coordinates
""" """
xmin, xpixel, _, ymax, _, ypixel = gdal_ds.GetGeoTransform() if isinstance(raster, pyotb.core.otbObject):
width, height = gdal_ds.RasterXSize, gdal_ds.RasterYSize info = pyotb.ReadImageInfo(raster)
xmax = xmin + width * xpixel spcx = info.GetParameterFloat('spacingx')
ymin = ymax + height * ypixel spcy = info.GetParameterFloat('spacingy')
orix = info.GetParameterFloat('originx')
oriy = info.GetParameterFloat('originy')
szx = info.GetParameterFloat('sizex')
szy = info.GetParameterFloat('sizey')
xmin = orix - .5 * spcx
ymax = oriy - .5 * spcy
else:
gdal_ds = gdal.Open(raster) if isinstance(raster, str) else raster
xmin, spcx, _, ymax, _, spcy = gdal_ds.GetGeoTransform()
szx, szy = gdal_ds.RasterXSize, gdal_ds.RasterYSize
xmax = xmin + szx * spcx
ymin = ymax + szy * spcy
return (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin) return (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)
def get_extent_wgs84(gdal_ds): def get_projection(raster):
""" """
Returns the extent in WGS84 CRS from a GDAL dataset Returns the projection (as str) of a raster.
Args: Args:
gdal_ds: GDAL dataset raster: GDAL dataset, str (raster filename) or otbobject
Returns:
a str
"""
if isinstance(raster, pyotb.core.otbObject):
return pyotb.ReadImageInfo(raster).GetParameterString("projectionref")
gdal_ds = gdal.Open(raster) if isinstance(raster, str) else raster
return gdal_ds.GetProjection()
def get_extent_wgs84(raster):
"""
Returns the extent in WGS84 CRS from a raster
Args:
raster: GDAL dataset, str (raster filename) or otbobject
Returns: Returns:
extent: coordinates in WGS84 CRS extent: coordinates in WGS84 CRS
""" """
coords = get_extent(gdal_ds) coords = get_extent(raster)
src_srs = osr.SpatialReference() src_srs = osr.SpatialReference()
src_srs.ImportFromWkt(gdal_ds.GetProjection()) projection = get_projection(raster)
src_srs.ImportFromWkt(projection)
tgt_srs = epsg2srs(4326) tgt_srs = epsg2srs(4326)
return reproject_coords(coords, src_srs, tgt_srs) return reproject_coords(coords, src_srs, tgt_srs)
...@@ -91,6 +122,5 @@ def get_bbox_wgs84(raster): ...@@ -91,6 +122,5 @@ def get_bbox_wgs84(raster):
a BoundingBox instance a BoundingBox instance
""" """
gdal_ds = gdal.Open(raster) if isinstance(raster, str) else raster extent = get_extent_wgs84(raster)
extent = get_extent_wgs84(gdal_ds)
return coord_list_to_bbox(extent) return coord_list_to_bbox(extent)
...@@ -3,14 +3,46 @@ from scenes_test_base import ScenesTestBase ...@@ -3,14 +3,46 @@ from scenes_test_base import ScenesTestBase
import tests_data import tests_data
from scenes import raster, vector from scenes import raster, vector
from scenes.spatial import BoundingBox from scenes.spatial import BoundingBox
import pyotb
import gdal
class SpatialTest(ScenesTestBase): class SpatialTest(ScenesTestBase):
BBOX_REF = BoundingBox(xmin=43.6279867026818,
xmax=43.70827724313909,
ymin=4.3153933040851165,
ymax=4.421982273366675)
def compare_bboxes(self, bbox1, bbox2, eps=1e-7):
coords = [(bbox1.xmin, bbox2.xmin),
(bbox1.xmax, bbox2.xmax),
(bbox1.ymin, bbox2.ymin),
(bbox1.ymax, bbox2.ymax)]
for v1, v2 in coords:
self.assertTrue(abs(v1 - v2) < eps)
def test_bbox(self): def test_bbox(self):
assert isinstance(raster.get_bbox_wgs84(tests_data.SCENE1.dimap_file_xs), BoundingBox) assert isinstance(raster.get_bbox_wgs84(tests_data.SCENE1.dimap_file_xs), BoundingBox)
assert isinstance(vector.get_bbox_wgs84(tests_data.ROI_MTP_4326), BoundingBox) assert isinstance(vector.get_bbox_wgs84(tests_data.ROI_MTP_4326), BoundingBox)
def test_otbinput_bbox(self):
img = pyotb.Input(tests_data.SCENE1.dimap_file_xs)
bbox_otb = raster.get_bbox_wgs84(img)
self.compare_bboxes(bbox_otb, self.BBOX_REF)
def test_otbapp_bbox(self):
img = pyotb.DynamicConvert(tests_data.SCENE1.dimap_file_xs)
bbox_otb = raster.get_bbox_wgs84(img)
self.compare_bboxes(bbox_otb, self.BBOX_REF)
def test_gdal_ds_bbox(self):
gdal_ds = gdal.Open(tests_data.SCENE1.dimap_file_xs)
bbox_gdal = raster.get_bbox_wgs84(gdal_ds)
self.compare_bboxes(bbox_gdal, self.BBOX_REF)
def test_filename_bbox(self):
bbox_file = raster.get_bbox_wgs84(tests_data.SCENE1.dimap_file_xs)
self.compare_bboxes(bbox_file, self.BBOX_REF)
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment