diff --git a/doc/index.md b/doc/index.md index cdc584c461daf79e714d6ebc152a23bee527b134..4b3d43eb5ef5a74f32ecfe7002d8d7f0ca521c3d 100644 --- a/doc/index.md +++ b/doc/index.md @@ -1,12 +1,14 @@ # Scenes -`Scenes` is a small python library aiming to provide a unified access to common remote sensing products. -Currently, supported sensors are: +`Scenes` is a small python library aiming to provide unified access to common +remote sensing products. +Currently, supported products are: -- Spot 6 -- Spot 7 +- Spot-6, Spot-7 (DIMAP) +- Spot-6, Spot-7 (IGN) - Sentinel-2 (Theia, Level 2) - Sentinel-2 (Theia, Level 3) +- Sentinel-2 (Microsoft Planetary Computer, Level 2) ## Installation @@ -16,25 +18,46 @@ To install the latest version: pip install https://gitlab.irstea.fr/umr-tetis/scenes/-/archive/develop/scenes-develop.zip --upgrade ``` -## Generic features +## Features + +### Sensor agnostic `Scene` instances have generic methods and attributes. - **Imagery sources**: imagery sources can be accessed using `get_xxx()` methods, where `xxx` is the name of the source to get. Any source can be used -in OTB, PyOTB pipelines, or with numpy directly. -- **Metadata**: metadata access with the `get_metadata()` method. +in OTB, PyOTB pipelines, or with numpy directly. Sources can be processed and +derived as new sources, the idea is to provide a set of methods for common +operations for a given sensor, and quickly chain processes from a particular +source. For instance, for Theia S2 Level-2 products, one can derive the +following `sc.get_b4().cld_msk(nodata=-10001)` consisting of the image +"drilled" with cloud mask. +- **Metadata**: metadata access with the `get_metadata()` method. It returns +a dict of information carried by the `Scene` instance. ``` py for key, value in sc.get_metadata(): print(f"md[{key}]: {value}") ``` -- **Acquisition date**: the `acquisition_date` member returns the acquisition date in the form of a `datetime.datetime` instance. -- **EPSG code**: the `epsg` member is an int for the EPSG of the scene. -- **Bounding box**: the `bbox_wgs84` member is a bounding box in WGS84 (EPSG 4326) coordinate reference system. +- **Acquisition date**: the `acquisition_date` member returns the acquisition +date in the form of a `datetime.datetime` instance. +- **EPSG code**: the `epsg` member is an `int` for the EPSG of the scene. +- **Bounding box**: the `bbox_wgs84` member is a bounding box in WGS84 (EPSG +4326) coordinate reference system. + +### STAC + +`Scene` instances can be created on the fly from `pystac.Item`. +See the STAC section in the API reference. + +## Supported sensors + +### Spot-6 and Spot-7 -## Spot-6 and Spot-7 +DIMAP products include cloud masks. This is not the case with IGN products. +IGN products cannot be processed with cloud drilling methods, and cannot be +optically calibrated. -### Instantiation +#### Instantiation Spot 6/7 scenes are instantiated from XS and PAN XML files of DIMAPS products. @@ -46,7 +69,7 @@ sc = Spot67Scene( ) ``` -### Image sources +#### Image sources Spot-6/7 scenes have 3 images sources: `xs`, `pan` and `pxs`. @@ -56,7 +79,7 @@ Spot-6/7 scenes have 3 images sources: `xs`, `pan` and `pxs`. | Pan | `get_pan()` | 1 channel (visible domain), ~1.5m physical spacing | | PXS | `get_pxs()` | 4 channels (same as XS), ~1.5m physical spacing | -Each source can be computed in DN/TOA/TOC reflectance. +Each source can be computed in DN/TOA/TOC reflectance: ``` py xs_raw = sc.get_xs() # DN @@ -77,10 +100,10 @@ The no-data value can be chosen. xs_toa_cld = xs_toa.cld_msk_drilled() # (1) ``` -1. To change the no-data inside the clouds mask: +To change the no-data inside the clouds mask: `masked_src = src.cld_msk_drilled(nodata=-999)` -### Examples +#### Examples Computing NDVI from XS image in TOA reflectance: @@ -128,13 +151,13 @@ superimposed = toa.resample_over(ref_img) superimposed.write("superimposed.tif") ``` -## Sentinel-2 +### Sentinel-2 Currently, Level-2 and Level-3 products from the [Theia land data center](https://www.theia-land.fr/en/product/sentinel-2-surface-reflectance/) and Level-2 products from Microsoft Planetary Computer are supported -### Instantiation +#### Instantiation Sentinel-2 scenes are instantiated from the product archive or folder. @@ -147,20 +170,17 @@ sc_level3 = Sentinel23AScene("/path/to/SENTINEL2X_...L3A_T31TEJ_D/.zip") Note that you can also instanciate uncomplete products, e.g. only 10m spacing bands. In this case, a warning will be displayed. -### Sources +#### Sources Sentinel-2 images include the following sources: + - single band, 10m spacing: `b4`, `b3`, `b2`, `b8` - single band, 20m spacing: `b5`, `b6`, `b7`, `b8a`, `b11`, `b12` - concatenated bands: `10m_bands`, `20m_bands` - for Theia products, quality masks: `clm_r1`, `clm_r2`, `edg_r1`, `edg_r2` -| Source | `Scene` method | Description | -|-----------|-------------------|-----------------------------------------------------------| -| 10m bands | `get_10m_bands()` | 10m physical spacing spectral bands (4, 3, 2, 8) | -| 20m bands | `get_20m_bands()` | 20m physical spacing spectral bands (5, 6, 7, 8a, 11, 12) | - -Each source can be drilled with no-data values. +For Theia products, imagery sources can be drilled from cloud mask (level 2) +or quality mask (level 3). For Level 2 products, using the cloud mask: @@ -169,7 +189,7 @@ src = sc_level2.get_10m_bands() masked_src = src.cld_msk_drilled() # (1) ``` -1. To change the no-data inside the clouds mask: +To change the no-data inside the clouds mask: `masked_src = src.cld_msk_drilled(nodata=-999)` For Level 3 products, using the quality mask: @@ -179,11 +199,35 @@ src = sc_level3.get_10m_bands() masked_src = src.flg_msk_drilled() # (1) ``` -1. To change the no-data inside the clouds mask: +To change the no-data inside the quality mask: `masked_src = src.flg_msk_drilled(nodata=-999)`. + Also, the `keep_flags_values` can be used to select the pixels to keep from a list of values in the quality mask (land, water, snow, ...). +## Bounding box, extent, epsg + +The bounding box in WGS84 coordinates of a particular raster, (py)otb pipeline +output, or source, can be computer using the following: + +```python +obj = ... +bbox = scenes.raster.get_bbox_wgs84(obj) +``` + +For the same kind of objects, the extent, in geo coordinates, can be retrieved +like this: + +```python +ext = scenes.raster.get_extent(obj) +``` + +The projection can also be retrieved: + +```python +proj = scenes.raster.get_projection(obj) +``` + ## Spatial and temporal indexation Scenes includes a module to perform spatial and temporal indexation of `Scene` @@ -224,7 +268,7 @@ matches3 = idx.find(bbox_wgs84=bbox, date_max="01/02/2019") # (3) keywords `date_min` and `date_max` -### Implement a new sensor from scratch +## Implement a new sensor from scratch You can implement quickly a new sensor in `scenes`. diff --git a/scenes/cache.py b/scenes/cache.py index b057f033e77dfd9d10467bdb765bb83ab8b054f4..a2fd2f00f514195b0f6aaf70d53da0cfee7e3ae5 100644 --- a/scenes/cache.py +++ b/scenes/cache.py @@ -38,10 +38,11 @@ class Cache(pyotb.Input): summary_modifier: optional function to modify the summary """ # Get app - pyotb_app = pyotb_output.pyotb_app + pyotb_app = pyotb_output.parent_pyotb_app # Get summary - summary = pyotb_app.summarize() # need pyotb >= 1.5.1 + assert isinstance(pyotb_app, pyotb.core.App) + summary = pyotb.summarize(pyotb_app) # need pyotb >= 1.5.1 # Modify summary if summary_modifier: @@ -56,7 +57,7 @@ class Cache(pyotb.Input): # Cache filename if not output_parameter_key: - output_parameter_key = pyotb_app.output_param + output_parameter_key = pyotb_app.output_image_key if not extension: extension = ".tif?&gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES" elif not extension.startswith("."): diff --git a/scenes/core.py b/scenes/core.py index df8b1be99a6a5dbdad466a03e59342922e427da7..42676b338cbfd174a8e75ea5382e221269807887 100644 --- a/scenes/core.py +++ b/scenes/core.py @@ -110,9 +110,9 @@ class Source(pyotb.Output): def __init__( self, root_scene: Scene, - out: Union[str, pyotb.core.otbObject], + out: Union[str, pyotb.core.OTBObject], parent: Source = None, - output_parameter_key: str = 'out' + output_parameter_key: str = None ): """ Args: @@ -120,11 +120,11 @@ class Source(pyotb.Output): out: image to deliver (can be an image filename (str), a pyotb.App, etc.) parent: parent Source instance - output_parameter_key: output parameter key of the app of `out` + output_parameter_key: output parameter key of the app """ - assert isinstance(root_scene, - Scene), f"root_scene type is {type(root_scene)}" + assert isinstance(root_scene, Scene), \ + f"root_scene type is {type(root_scene)}" self.root_scene = root_scene # Here we call the pyotb.Output() constructor. # Since it can only be called with pyotb apps, we do the following: @@ -135,9 +135,12 @@ class Source(pyotb.Output): if isinstance(out, str): app = pyotb.Input(out) elif isinstance(out, pyotb.Output): - app = out.pyotb_app + app = out.parent_pyotb_app - super().__init__(app=app, output_parameter_key=output_parameter_key) + super().__init__( + pyotb_app=app, + param_key=output_parameter_key or app.output_image_key + ) assert parent is not self, "You cannot assign a new source to its " \ "parent instance" self.parent = parent # parent source (is another Source instance) @@ -159,15 +162,18 @@ class Source(pyotb.Output): """ for new_app in args: self._app_stack.append(new_app) - return self.__class__(root_scene=self.root_scene, - out=self._app_stack[-1], parent=self, **kwargs) + return self.__class__( + root_scene=self.root_scene, + out=self._app_stack[-1], + parent=self, **kwargs + ) def __getattr__(self, name: str): """ This function is called when an attribute or a method has been called, but not found in the object properties. We override it to avoid falling back into the depths of - `pyotb.otbObject`. + `pyotb.OTBObject`. Args: name: name of the attribute or method to access @@ -187,7 +193,7 @@ class Source(pyotb.Output): a string summarizing the metadata """ - return pprint(self.summarize()) + return pprint(pyotb.summarize(self.parent_pyotb_app)) class CommonImagerySource(Source): @@ -234,7 +240,7 @@ class CommonImagerySource(Source): def masked( self, - binary_mask: Union[str, pyotb.core.otbObject], + binary_mask: Union[str, pyotb.core.OTBObject], nodata: Union[float, int] = 0 ) -> CommonImagerySource: """ @@ -260,7 +266,7 @@ class CommonImagerySource(Source): def resample_over( self, - ref_img: Union[str, pyotb.core.otbObject], + ref_img: Union[str, pyotb.core.OTBObject], interpolator: str = "bco", nodata: Union[float, int] = 0 ) -> CommonImagerySource: @@ -310,7 +316,7 @@ class CommonImagerySource(Source): def clip_over_img( self, - ref_img: Union[str, pyotb.core.otbObject] + ref_img: Union[str, pyotb.core.OTBObject] ) -> CommonImagerySource: """ Return the source clipped over the ROI specified by the input image diff --git a/scenes/deepnets.py b/scenes/deepnets.py index 656cf1a2da3066cdb7e5fac030ecb4085191b281..7a88e6016c6b7977e52099d90f8649c03ce72f69 100644 --- a/scenes/deepnets.py +++ b/scenes/deepnets.py @@ -18,7 +18,7 @@ import scenes archive = "SENTINEL2B_..._T31TEJ_D_V1-8.zip" s2_scene = scenes.sentinel.Sentinel22AScene(archive) s2_image = s2_scene.get_10m_bands() -sr = scenes.deepnets.sr4rs(s2_image) # pyotb.core.otbObject +sr = scenes.deepnets.sr4rs(s2_image) # pyotb.core.OTBObject sr.write("sr.tif") ``` """ @@ -56,7 +56,7 @@ def inference(dic: Dict[str, any]): def sr4rs( - input_image: Union[str, pyotb.core.otbObject], + input_image: Union[str, pyotb.core.OTBObject], model_url: str = SR4RS_MODEL_URL, tmp_dir: str = "/tmp" ) -> pyotb.core.App: diff --git a/scenes/raster.py b/scenes/raster.py index af09efb2c26b2851d3b00dff9992a7ff0f08a0e9..ff94b4681a56b31eae1072ae6d3a6c2a17054622 100644 --- a/scenes/raster.py +++ b/scenes/raster.py @@ -28,26 +28,26 @@ def get_epsg(gdal_ds: gdal.Dataset) -> int: def get_extent( - raster: Union[gdal.Dataset, str, pyotb.core.otbObject] + raster: Union[gdal.Dataset, str, pyotb.core.OTBObject] ) -> Tuple[Tuple[float]]: """ Return list of corners coordinates from a raster Args: - raster: GDAL dataset, str (raster filename) or otbobject + raster: GDAL dataset, str (raster filename) or OTBObject Returns: list of coordinates """ - if isinstance(raster, pyotb.core.otbObject): + 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') + spcx = info['spacingx'] + spcy = info['spacingy'] + orix = info['originx'] + oriy = info['originy'] + szx = info['sizex'] + szy = info['sizey'] xmin = orix - .5 * spcx ymax = oriy - .5 * spcy else: @@ -61,32 +61,32 @@ def get_extent( def get_projection( - raster: Union[gdal.Dataset, str, pyotb.core.otbObject] + raster: Union[gdal.Dataset, str, pyotb.core.OTBObject] ) -> str: """ Returns the projection (as str) of a raster. Args: - raster: GDAL dataset, str (raster filename) or otbobject + raster: GDAL dataset, str (raster filename) or OTBObject Returns: a str """ - if isinstance(raster, pyotb.core.otbObject): - return pyotb.ReadImageInfo(raster).GetParameterString("projectionref") + if isinstance(raster, pyotb.core.OTBObject): + return pyotb.ReadImageInfo(raster)["projectionref"] gdal_ds = gdal.Open(raster) if isinstance(raster, str) else raster return gdal_ds.GetProjection() def get_extent_wgs84( - raster: Union[gdal.Dataset, str, pyotb.core.otbObject] + raster: Union[gdal.Dataset, str, pyotb.core.OTBObject] ) -> Tuple[Tuple[float]]: """ Returns the extent in WGS84 CRS from a raster Args: - raster: GDAL dataset, str (raster filename) or otbobject + raster: GDAL dataset, str (raster filename) or OTBObject Returns: extent: coordinates in WGS84 CRS @@ -121,7 +121,7 @@ def get_epsg_extent_wgs84(filename: str) -> Tuple[int, Tuple[Tuple[float]]]: def get_bbox_wgs84( - raster: Union[gdal.Dataset, str, pyotb.core.otbObject] + raster: Union[gdal.Dataset, str, pyotb.core.OTBObject] ) -> BoundingBox: """ Returns the bounding box from a raster. diff --git a/scenes/sentinel.py b/scenes/sentinel.py index 469a57abc13e6f9e0ee335f312757c185ac0355f..792efafd3a33f7178d6a9e5aad692549b03c6699 100644 --- a/scenes/sentinel.py +++ b/scenes/sentinel.py @@ -174,7 +174,7 @@ Note that the resulting transformed `Sentinel2Theia2ASource` and ## Usage with pyotb -As `scenes.core.Source`, it also can be used like any `pyotb.core.otbObject`. +As `scenes.core.Source`, it also can be used like any `pyotb.core.OTBObject`. The following example show how to use an OTB application with a source at input. @@ -218,7 +218,7 @@ class Sentinel2TheiaGenericSource(CommonImagerySource): new masked source """ - img_spc = pyotb.ReadImageInfo(self).GetParameterInt('spacingx') + img_spc = pyotb.ReadImageInfo(self)['spacingx'] if img_spc not in msk_dict: raise KeyError(f"No mask for image spacing {img_spc}") binary_mask = pyotb.BandMath({"il": msk_dict[img_spc], "exp": exp}) diff --git a/scenes/spot.py b/scenes/spot.py index 6d2792a819db273127d3c08588ca181c68e7095e..0f97e39f4e342e170c3ce02c35cbf115f3eb1234 100644 --- a/scenes/spot.py +++ b/scenes/spot.py @@ -128,7 +128,7 @@ Note that the resulting transformed `Spot67DRSSource` is still an instance of # Usage with pyotb -As `scenes.core.Source`, it also can be used like any `pyotb.core.otbObject`. +As `scenes.core.Source`, it also can be used like any `pyotb.core.OTBObject`. The following example show how to use an OTB application with a source at input. diff --git a/setup.py b/setup.py index 886285c35ebca844602d434d290dc808efb8b9dc..3407b82bdde6988d4f7e93355e5d5a7d2757eb6e 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ with open("README.md", "r", encoding="utf-8") as fh: setuptools.setup( name="scenes", - version="1.6.4", + version="1.7.0", author="Rémi Cresson", author_email="remi.cresson@inrae.fr", description="Library to ease the processing of remote sensing products", @@ -17,7 +17,7 @@ setuptools.setup( python_requires=">=3.6", install_requires=[ "rtree", - "pyotb==1.5.4", + "pyotb==2.0.0-dev4", "requests", "tqdm", "pystac-client", diff --git a/test/scenes_test_base.py b/test/scenes_test_base.py index 1e0be36527c2ae8191ff8e80843d54898a06ba30..fa9b83e8ab18524348723248b8c2b26b27500f22 100644 --- a/test/scenes_test_base.py +++ b/test/scenes_test_base.py @@ -1,10 +1,9 @@ # -*- coding: utf-8 -*- import filecmp -import pyotb import unittest from abc import ABC -from scenes import utils +import pyotb class ScenesTestBase(ABC, unittest.TestCase): @@ -21,10 +20,15 @@ class ScenesTestBase(ABC, unittest.TestCase): :param mae_threshold: mean absolute error threshold :return: boolean """ - nbands_in = pyotb.Input(image).shape[-1] - nbands_ref = pyotb.Input(reference).shape[-1] + if isinstance(image, str): + image = pyotb.Input(image) + if isinstance(reference, str): + reference = pyotb.Input(reference) + nbands_in = image.shape[-1] + nbands_ref = reference.shape[-1] - assert nbands_in == nbands_ref, f"Image has {nbands_in} but baseline has {nbands_ref}" + assert nbands_in == nbands_ref, \ + f"Image has {nbands_in} but baseline has {nbands_ref}" for i in range(1, nbands_ref + 1): dic = { @@ -41,10 +45,8 @@ class ScenesTestBase(ABC, unittest.TestCase): "roi.sizey": roi[3] }) comp = pyotb.CompareImages(dic) - mae = comp.GetParameterFloat('mae') - assert mae < mae_threshold, f"MAE is {mae} > {mae_threshold}\n" \ - f"\nTarget summary:\n{utils.pprint(image.summarize())}\n" \ - f"\nReference summary:\n{utils.pprint(reference.summarize())}" + mae = comp['mae'] + assert mae < mae_threshold, f"MAE is {mae} > {mae_threshold}\n" def compare_file(self, file, reference): """ @@ -57,5 +59,7 @@ class ScenesTestBase(ABC, unittest.TestCase): Return: a boolean """ - assert filecmp.cmp(file, - reference), f"File {file} differs with baseline ({reference})" + assert filecmp.cmp( + file, + reference + ), f"File {file} differs with baseline ({reference})" diff --git a/test/spot67_imagery_test.py b/test/spot67_imagery_test.py index 6d1afbf7c5181f5b20be53de901a3fac9dfe0213..d625833266c5f661e41ac63b78505452d8df2b85 100644 --- a/test/spot67_imagery_test.py +++ b/test/spot67_imagery_test.py @@ -31,7 +31,7 @@ pxs = sc.get_pxs() @singledispatch def subset(inp: Any) -> Any: raise TypeError( - "Invalid type, must be one of: pyotb.otbObject, Source" + "Invalid type, must be one of: pyotb.core.OTBObject, Source" ) @@ -46,7 +46,7 @@ def subset_source(inp: scenes.core.Source): @subset.register -def subset_pyotb(inp: pyotb.otbObject): +def subset_pyotb(inp: pyotb.core.OTBObject): return inp[roi[0]:roi[0] + roi[2], roi[1]:roi[1] + roi[3], :]