diff --git a/apps/drs_spot67_import.py b/apps/drs_spot67_import.py index b0c9d58eed2fb03202a2937fdbee2c92257e3eac..d8b5ec59149673729461cf725a81016f6b907636 100755 --- a/apps/drs_spot67_import.py +++ b/apps/drs_spot67_import.py @@ -6,7 +6,6 @@ from scenes.spot import get_spot67_scenes # Arguments 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("--tile_name", help="Filter only images from this tile (optional)") parser.add_argument("--out_pickle", help="Output pickle file", required=True) params = parser.parse_args() diff --git a/scenes/raster.py b/scenes/raster.py index 0c4c1685d497bd16ba86525a5e383120e572fec9..6cf66bd55c5926ef65ae312d77a2b3a749b68eaa 100644 --- a/scenes/raster.py +++ b/scenes/raster.py @@ -1,6 +1,7 @@ """ This module contains a set of functions to deal with GDAL datasets (rasters) """ +import pyotb from osgeo import osr, gdal from scenes.vector import reproject_coords from scenes.spatial import epsg2srs, coord_list_to_bbox @@ -23,39 +24,69 @@ def get_epsg(gdal_ds): 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: - gdal_ds: GDAL dataset + raster: GDAL dataset, str (raster filename) or otbobject Returns: list of coordinates """ - xmin, xpixel, _, ymax, _, ypixel = gdal_ds.GetGeoTransform() - width, height = gdal_ds.RasterXSize, gdal_ds.RasterYSize - xmax = xmin + width * xpixel - ymin = ymax + height * ypixel + if isinstance(raster, pyotb.core.otbObject): + info = pyotb.ReadImageInfo(raster) + spcx = info.GetParameterFloat('spacingx') + 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) -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: - 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: extent: coordinates in WGS84 CRS """ - coords = get_extent(gdal_ds) + coords = get_extent(raster) src_srs = osr.SpatialReference() - src_srs.ImportFromWkt(gdal_ds.GetProjection()) + projection = get_projection(raster) + src_srs.ImportFromWkt(projection) tgt_srs = epsg2srs(4326) return reproject_coords(coords, src_srs, tgt_srs) @@ -91,6 +122,5 @@ def get_bbox_wgs84(raster): a BoundingBox instance """ - gdal_ds = gdal.Open(raster) if isinstance(raster, str) else raster - extent = get_extent_wgs84(gdal_ds) + extent = get_extent_wgs84(raster) return coord_list_to_bbox(extent) diff --git a/test/spatial_test.py b/test/spatial_test.py index 99247be4b7dd02b2432c2c79da646f9ba108cf57..a05bc5fba078f53b4a10db8066895a0ec6252117 100644 --- a/test/spatial_test.py +++ b/test/spatial_test.py @@ -3,14 +3,46 @@ from scenes_test_base import ScenesTestBase import tests_data from scenes import raster, vector from scenes.spatial import BoundingBox +import pyotb +import gdal 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): 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) + 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__': unittest.main()