Scenes¶

Scenes is a small python library aiming to provide unified access to common remote sensing products. Currently, supported products are:

• 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¶

To install the latest version:

pip install https://gitlab.irstea.fr/umr-tetis/scenes/-/archive/develop/scenes-develop.zip --upgrade

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. 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.

  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.

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¶

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¶

Spot 6/7 scenes are instantiated from XS and PAN XML files of DIMAPS products.

from scenes.spot import Spot67Scene
sc = Spot67Scene(
    dimap_file_xs="/path/to/DIM_SPOT7_MS_..._1.XML",
    dimap_file_pan="/path/to/DIM_SPOT7_P_..._1.XML"
)

Image sources¶

Spot-6/7 scenes have 3 images sources: xs, pan and pxs.

Each source can be computed in DN/TOA/TOC reflectance:

xs_raw = sc.get_xs()  # DN
xs_toa = sc.get_xs().reflectance(level="toa")  # TOA
xs_toc = sc.get_xs().reflectance(level="toc")  # TOC

!!! Warning

OTB version must be >= 8.1.0 to support optical calibration using source 
`reflectance()` method.
kwargs of `reflectance()` are the same as `pyotb.OpticalCalibration`.

Each source can be masked with the cloud masks of the original product. The no-data value can be chosen.

xs_toa_cld = xs_toa.cld_msk_drilled()  # (1)

To change the no-data inside the clouds mask: masked_src = src.cld_msk_drilled(nodata=-999)

Examples¶

Computing NDVI from XS image in TOA reflectance:

xs = toa.get_xs().reflectance(level="toa")  # (1)
exp = "(im1b4-im1b1)/(im1b4+im1b1)"
ndvi = pyotb.bandmath(exp=exp, il=[xs])     # (2)
ndvi.write("ndvi.tif")

1. xs is a scenes.spot.Spot67Source instance
2. ndvi is a pyotb.App object that inputs xs

The next example is a set of preprocessing operations on a Spot-6/7 XS image:

1. TOA reflectance
2. pansharpening
3. replacing the pixels under the clouds masks with "0"
4. clip the result over a reference raster

pxs = sc.get_pxs().reflectance(level="toa")  # (1)
drilled = pxs.cld_msk_drilled()              # (2)
ref_img = "/tmp/S2A_2020...._FRE_10m.tif"
subset = drilled.clip_over_img(ref_img)      # (3)
subset.write("subset.tif")

1. pxs is a scenes.spot.Spot67Source instance
2. drilled is a scenes.spot.Spot67Source instance
3. subset is a scenes.spot.Spot67Source instance

Note that we could have changed the no-data value in the cloud masking:

drilled = pxs.drilled(nodata=-999)

Superimpose an image over a reference image. In the example below, ref_img is another scenes.core.Source instance.

toa = sc.get_pxs().reflectance(level="toa")
superimposed = toa.resample_over(ref_img)
superimposed.write("superimposed.tif")

Sentinel-2¶

Currently, Level-2 and Level-3 products from the Theia land data center and Level-2 products from Microsoft Planetary Computer are supported

Instantiation¶

Sentinel-2 scenes are instantiated from the product archive or folder.

from scenes.sentinel import Sentinel22AScene, Sentinel23AScene
sc_level2 = Sentinel22AScene("/path/to/SENTINEL2A_..._V1-8/.zip")
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¶

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

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:

src = sc_level2.get_10m_bands()
masked_src = src.cld_msk_drilled()  # (1)

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:

src = sc_level3.get_10m_bands()
masked_src = src.flg_msk_drilled()  # (1)

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:

obj = ...
bbox = scenes.raster.get_bbox_wgs84(obj)

For the same kind of objects, the extent, in geo coordinates, can be retrieved like this:

ext = scenes.raster.get_extent(obj)

The projection can also be retrieved:

proj = scenes.raster.get_projection(obj)

Spatial and temporal indexation¶

Scenes includes a module to perform spatial and temporal indexation of Scene instances.

Query in space¶

Perform a query in space (WGS84 bounding box) and time (optional) with an indexation structure.

from scenes import spot, indexation
scs = spot.get_spot67_scenes(root_dir)  # (1)
idx = indexation.Index(scenes)          # (2)
matches = idx.find(bbox_wgs84=bbox)     # (3)
for sc in matches:
  print(sc)

1. scs is a list of scenes.core.Scene instances
2. idx is a scenes.indexation.Index instance, namely a spatio-temporal index
3. matches is a list of scenes.core.Scene instances

Query in space and time¶

To perform searches in time:

matches1 = idx.find(bbox_wgs84=bbox, "01/02/2019", "01-12-2020")  # (1)
matches2 = idx.find(bbox_wgs84=bbox, "01/02/2019")                # (2)
matches3 = idx.find(bbox_wgs84=bbox, date_max="01/02/2019")       # (3)

1. A few date formats are supported: "dd/mm/yyyy", "dd-mm-yyyy", and "yyyy-mm-dd"
2. Here we don't have upper bound
3. You can also specify only an upper bound, without lower bound (use the keywords date_min and date_max

Implement a new sensor from scratch¶

You can implement quickly a new sensor in scenes.

• One or multiple MySensorSource1, ..., MySensorSourceM classes, inheriting from scenes.core.Source, and implementing the image access for the new sensor.
• A new MySensorScene class, inheriting from scenes.core.Scene. This class must provide one or multiple methods to its sources.