Mosaic

Basic example of a mosaic.

Initialization of the notebook

• Imports

• Loggers

• Paths

import os

from cloudpathlib import AnyPath
from eosets.mosaic import Mosaic
from eoreader.bands import MNDWI, RED, SLOPE, CLOUDS
from eoreader.env_vars import DEM_PATH

# Create logger
import logging
from sertit import logs

logs.init_logger(logging.getLogger("eoreader"))
logs.init_logger(logging.getLogger("eosets"))

# Get the base paths
data_path = AnyPath(r"/home/ds2_db3/CI/eosets/DATA")
db_path = AnyPath(r"/home/ds2_db2/BASES_DE_DONNEES")

# Set a DEM in order to compte the slope over the mosaic
dem_sub_dir_path = ["GLOBAL", "COPDEM_30m", "COPDEM_30m.vrt"]
os.environ[DEM_PATH] = str(db_path.joinpath(*dem_sub_dir_path))

# Get a list of Sentinel-2 data with the same day
s2_paths = [
    data_path / "S2B_MSIL2A_20220228T102849_N0400_R108_T32TLT_20220228T134712.SAFE",
    data_path / "S2B_MSIL2A_20220228T102849_N0400_R108_T32ULU_20220228T134712.SAFE",
    data_path / "S2B_MSIL2A_20220228T102849_N0400_R108_T32ULV_20220228T134712.SAFE",
]

Creation of the mosaic

You just have to give the paths to your Sentinel-2 producs. Their extent or their footprints should be contiguous, and they need to be acquired on the same day.

You can choose a mosaicing method. VRT advised as it is lightweight and faster.

# Then with compatible
mosaic = Mosaic(s2_paths, mosaic_method="VRT")

Stack

Create a stack of 4 bands, mosaicing the 3 Sentinel-2 products.

# Stack with a resolution of 60m
stack = mosaic.stack(
    [MNDWI, RED, SLOPE, CLOUDS],
    pixel_size=60,
)

2024-04-25 09:00:21,592 - [DEBUG] - *** Loading ['MNDWI', 'RED', 'SLOPE', 'CLOUDS'] for 20220228T102849_S2_T32TLT_L2A_134712 ***

2024-04-25 09:00:21,625 - [DEBUG] - Loading bands ['RED', 'GREEN', 'SWIR_1']

2024-04-25 09:00:21,861 - [DEBUG] - Read RED

2024-04-25 09:00:22,204 - [DEBUG] - Manage nodata for band RED

2024-04-25 09:00:23,124 - [DEBUG] - Converting RED to reflectance

2024-04-25 09:00:34,236 - [DEBUG] - Read GREEN

2024-04-25 09:00:34,524 - [DEBUG] - Manage nodata for band GREEN

2024-04-25 09:00:34,969 - [DEBUG] - Converting GREEN to reflectance

2024-04-25 09:00:39,494 - [DEBUG] - Read SWIR_1

2024-04-25 09:00:39,661 - [DEBUG] - Manage nodata for band SWIR_1

2024-04-25 09:00:39,875 - [DEBUG] - Converting SWIR_1 to reflectance

2024-04-25 09:00:42,202 - [DEBUG] - Loading indices ['MNDWI']

2024-04-25 09:00:42,266 - [DEBUG] - Loading DEM bands ['SLOPE']

2024-04-25 09:00:42,266 - [DEBUG] - Warping DEM for 20220228T102849_S2_T32TLT_L2A_134712

2024-04-25 09:00:42,285 - [DEBUG] - Using DEM: /home/ds2_db2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt

2024-04-25 09:00:44,908 - [DEBUG] - Computing slope for 20220228T102849_S2_T32TLT_L2A_134712

2024-04-25 09:01:18,194 - [DEBUG] - Loading Cloud bands ['CLOUDS']

2024-04-25 09:01:19,271 - [DEBUG] - *** Loading ['MNDWI', 'RED', 'SLOPE', 'CLOUDS'] for 20220228T102849_S2_T32ULU_L2A_134712 ***

2024-04-25 09:01:19,287 - [DEBUG] - Loading bands ['RED', 'GREEN', 'SWIR_1']

2024-04-25 09:01:19,455 - [DEBUG] - Read RED

2024-04-25 09:01:19,631 - [DEBUG] - Manage nodata for band RED

2024-04-25 09:01:20,082 - [DEBUG] - Converting RED to reflectance

2024-04-25 09:01:25,013 - [DEBUG] - Read GREEN

2024-04-25 09:01:25,104 - [DEBUG] - Manage nodata for band GREEN

2024-04-25 09:01:25,400 - [DEBUG] - Converting GREEN to reflectance

2024-04-25 09:01:29,572 - [DEBUG] - Read SWIR_1

2024-04-25 09:01:29,783 - [DEBUG] - Manage nodata for band SWIR_1

2024-04-25 09:01:30,052 - [DEBUG] - Converting SWIR_1 to reflectance

2024-04-25 09:01:32,821 - [DEBUG] - Loading indices ['MNDWI']

2024-04-25 09:01:32,867 - [DEBUG] - Loading DEM bands ['SLOPE']

2024-04-25 09:01:32,868 - [DEBUG] - Warping DEM for 20220228T102849_S2_T32ULU_L2A_134712

2024-04-25 09:01:32,887 - [DEBUG] - Using DEM: /home/ds2_db2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt

2024-04-25 09:01:35,220 - [DEBUG] - Computing slope for 20220228T102849_S2_T32ULU_L2A_134712

2024-04-25 09:01:46,217 - [DEBUG] - Loading Cloud bands ['CLOUDS']

2024-04-25 09:01:47,286 - [DEBUG] - *** Loading ['MNDWI', 'RED', 'SLOPE', 'CLOUDS'] for 20220228T102849_S2_T32ULV_L2A_134712 ***

2024-04-25 09:01:47,306 - [DEBUG] - Loading bands ['RED', 'GREEN', 'SWIR_1']

2024-04-25 09:01:47,544 - [DEBUG] - Read RED

2024-04-25 09:01:47,792 - [DEBUG] - Manage nodata for band RED

2024-04-25 09:01:48,300 - [DEBUG] - Converting RED to reflectance

2024-04-25 09:01:52,840 - [DEBUG] - Read GREEN

2024-04-25 09:01:52,988 - [DEBUG] - Manage nodata for band GREEN

2024-04-25 09:01:53,220 - [DEBUG] - Converting GREEN to reflectance

2024-04-25 09:01:57,377 - [DEBUG] - Read SWIR_1

2024-04-25 09:01:57,513 - [DEBUG] - Manage nodata for band SWIR_1

2024-04-25 09:01:57,800 - [DEBUG] - Converting SWIR_1 to reflectance

2024-04-25 09:02:00,187 - [DEBUG] - Loading indices ['MNDWI']

2024-04-25 09:02:00,233 - [DEBUG] - Loading DEM bands ['SLOPE']

2024-04-25 09:02:00,234 - [DEBUG] - Warping DEM for 20220228T102849_S2_T32ULV_L2A_134712

2024-04-25 09:02:00,251 - [DEBUG] - Using DEM: /home/ds2_db2/BASES_DE_DONNEES/GLOBAL/COPDEM_30m/COPDEM_30m.vrt

2024-04-25 09:02:02,690 - [DEBUG] - Computing slope for 20220228T102849_S2_T32ULV_L2A_134712

2024-04-25 09:02:12,779 - [DEBUG] - Loading Cloud bands ['CLOUDS']

2024-04-25 09:02:13,880 - [DEBUG] - Merging bands MNDWI

2024-04-25 09:02:14,009 - [DEBUG] - Merging bands RED

2024-04-25 09:02:14,137 - [DEBUG] - Merging bands SLOPE

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 09:02:14,267 - [DEBUG] - Merging bands CLOUDS

2024-04-25 09:02:14,394 - [DEBUG] - Collocating bands

0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 09:02:15,873 - [DEBUG] - Stacking

/opt/conda/lib/python3.9/site-packages/xarray/core/indexes.py:662: RuntimeWarning: '<' not supported between instances of 'DemBandNames' and 'SpectralBandNames', sort order is undefined for incomparable objects.
  new_pd_index = pd_indexes[0].append(pd_indexes[1:])

# Scale the mosaic in order to print it
from sertit.display import scale
scaled = stack.copy(data=scale(stack.data))

/opt/conda/lib/python3.9/site-packages/dask/array/core.py:1712: FutureWarning: The `numpy.nanpercentile` function is not implemented by Dask array. You may want to use the da.map_blocks function or something similar to silence this warning. Your code may stop working in a future release.
  warnings.warn(

# Plot the mosaic
import cartopy.crs as ccrs
crs = ccrs.UTM("32")
scaled[0:3, :, :].plot.imshow(
    robust=True,
    transform=crs,
    x="x",
    y="y",
    figsize=[10, 30],
    subplot_kws={'projection':crs}
)

<matplotlib.image.AxesImage at 0x7f5c323ffee0>