9. Sentinel 2 Mosaics¶

This tutorial covers the complete pipeline for producing seamless mosaics of Sentinel 2 imagery using the s2mosaic module from eumap and AWS.

This is done in 3 phases: 1. For each S2 tile that overlaps the ROI, a composite is made from all images that fall in designated time range. This composite can be the median or any other percentile (25%, 75%, …), but all pixels that are marked as cloud or cloud shadow are excluded first. 2. A mosaic of all tiles that are acquired along the same orbital track is made (overlapping pixels from different tiles are averaged), for each orbital track that intersects the ROI. 3. The orbit-wise mosaics are stitched together by calculating a weighted mean along regions where they overlap. The distance of each pixel from their respective orbital track is taken as the weight.

For this tutorial to work you will need to register for an AWS account and configure AWS CLI tools: https://docs.aws.amazon.com/cli/latest/userguide/cli-chap-configure.html

[ ]:

import eumap
import eumap.datasets.eo.s2mosaic
from eumap.datasets.eo.s2mosaic.satmos_client import JobSchedulerLocal, SatMosClient
from eumap.datasets.eo.s2mosaic.mosaic_maker import MosaicMaker

from datetime import datetime

9.1. 1. S2 tile composites¶

To complete the compositing phase we need to define a time period and obtain a list of S2 scenes as well as a list of the S2 tiles that overlap our ROI. Here we will set some temporary folders for intermediate results and define a timeframe and spatial coverage of mosaic.

[5]:

source = 's2l2a'
name = 'summer2019'
band = 'B08'                                        # we will calculate only the blue S2 band

scenes_csv = '09_scenes_tutorial.csv'               # file with all scenes (images)
folder_data = '/data/gh/tutorial/data'              # temporary folder for downloading scenes
folder_tmp = '/data/gh/tutorial/tmp'                # temporary folder for creating tiffs
folder_tiles = '/data/gh/tutorial/tiles'            # output folder for tiles composites
bucket = '/data/gh/'                                # name of the output s3 bucket
folder_out_parent = 'tutorial/tiles'

mosaic_name = f'{source}_{name}_{band}'             # name of the mosaic
date_from = datetime(2019, 6, 21)
date_to = datetime(2019, 9, 23)
s2tiles = ['32UME','32UNE','32UPE','32UQE','33UUV', # these S2 tiles define area of mosaic
           '32UMD','32UND','32UPD','32UQD','33UUU',
           '32UMC','32UNC','32UPC','32UQC','33UUT']
nworkers = 4                                        # number of parallel workers

The JobSchedulerLocal object takes care of queueing jobs, while SatMosClient takes care of parallelization.

[7]:

job_scheduler = JobSchedulerLocal(
    source=source,
    out_folder_prefix = mosaic_name,
    band = band,
    scenes_csv = scenes_csv,
    from_date = date_from, to_date=date_to,
    bucket = bucket,
    debug = True,
    data_folder = folder_data,
    tmp_folder = folder_tmp,
    out_parent_folder = folder_out_parent,
    tiles = s2tiles
)

client = SatMosClient(job_scheduler, nworkers=nworkers)

Once we have all of this defined we can run the client. For the job defined above the processing can take up to 3 hours. Results are tiles in /data/gh/tutorial/s2l2a_summer2019_B08. There is 4 folders: CNT, P25, P50, P75. CNT means number of source images for each pixel (excluding clouds), while the rest are the 25th, 50th and 75th percentiles.

[8]:

client.run()

STARTING: index=1, s2l2a_summer2019_B08, 32UME
New job: avail:58.1074333190918 GB, n_jobs: 1
STARTING: index=2, s2l2a_summer2019_B08, 32UNE
New job: avail:54.0652961730957 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=3, s2l2a_summer2019_B08, 32UPE
New job: avail:41.1991081237793 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=4, s2l2a_summer2019_B08, 32UQE
New job: avail:57.94245529174805 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=5, s2l2a_summer2019_B08, 33UUV
New job: avail:46.87162780761719 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=6, s2l2a_summer2019_B08, 32UMD
New job: avail:58.21084976196289 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=7, s2l2a_summer2019_B08, 32UND
New job: avail:52.733768463134766 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=8, s2l2a_summer2019_B08, 32UPD
New job: avail:58.11817169189453 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=9, s2l2a_summer2019_B08, 32UQD
New job: avail:46.91321563720703 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=10, s2l2a_summer2019_B08, 33UUU
New job: avail:58.18669509887695 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=11, s2l2a_summer2019_B08, 32UMC
New job: avail:47.37238311767578 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=12, s2l2a_summer2019_B08, 32UNC
New job: avail:58.22690963745117 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=13, s2l2a_summer2019_B08, 32UPC
New job: avail:47.235713958740234 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=14, s2l2a_summer2019_B08, 32UQC
New job: avail:58.154666900634766 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
STARTING: index=15, s2l2a_summer2019_B08, 33UUT
New job: avail:47.84326934814453 GB, n_jobs: 2
FINISHED TEST: success=True, report=OK
FINISHED TEST: success=True, report=OK
No jobs for me !!!

9.1.1. 2. Relative orbits¶

S2 scenes that are acquired along the same orbital track are very likely imaged in similar conditions so we can mosaic them by simply averaging the overlaps.

[12]:

bucket = None   # input tiles are stored locally
folder_tmp = '/data/gh/tutorial/tmp'

nworkers = 20

dxy = 30        #output resolution
crs = '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs'
resampling_method = 'average'   # 'average', 'nearest', 'cubic'

[13]:

mosaic_maker = MosaicMaker(
    nworkers=nworkers,
    dst_crs = crs,
    dxy = dxy,
    bucket = bucket,
    tmp_folder=folder_tmp,
    debug=True
)

Now we can run MosaicMaker, which takes composites from the tiles folder and stitches together ones along the same orbital track. This may take a few minutes to finish, depending on ‘nworkers’.

[17]:

for p in ['CNT','P25','P50','P75']:
    folder_input = f'/data/gh/tutorial/tiles/s2l2a_summer2019_B08/{p}/'
    folder_output = f'/data/gh/tutorial/orbits/s2l2a_summer2019_B08/{p}/'
    print(folder_input)
    mosaic_maker.mosaic_orbits_tiled(folder_input, folder_output, resampling_method=resampling_method)

/data/gh/tutorial/tiles/s2l2a_summer2019_B08/CNT/
[2022-06-29 18:10:40] mosaic: /data/gh/tutorial/tiles/s2l2a_summer2019_B08/CNT/
[2022-06-29 18:10:40] number of orbits: 5

1/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R051.tif, 1 files:  20%|██        | 1/5 [00:00<00:00, 5890.88it/s]


Window 0: 100%|██████████| 1/1 [00:01<00:00,  1.89s/it]
2/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R022.tif, 6 files:  40%|████      | 2/5 [00:02<00:03,  1.00s/it]




Window 1: 100%|██████████| 2/2 [00:07<00:00,  3.82s/it]
3/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R008.tif, 6 files:  60%|██████    | 3/5 [00:10<00:07,  3.97s/it]




Window 1: 100%|██████████| 2/2 [00:07<00:00,  3.82s/it]
4/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R065.tif, 12 files:  80%|████████  | 4/5 [00:18<00:05,  5.50s/it]








Window 1: 100%|██████████| 4/4 [00:11<00:00,  2.99s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R108.tif, 13 files: 100%|██████████| 5/5 [00:31<00:00,  8.06s/it]








Window 1: 100%|██████████| 4/4 [00:15<00:00,  3.87s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R108.tif, 13 files: 100%|██████████| 5/5 [00:47<00:00,  9.51s/it]

/data/gh/tutorial/tiles/s2l2a_summer2019_B08/P25/
[2022-06-29 18:11:28] mosaic: /data/gh/tutorial/tiles/s2l2a_summer2019_B08/P25/
[2022-06-29 18:11:28] number of orbits: 5


1/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R051.tif, 1 files:  20%|██        | 1/5 [00:00<00:00, 10058.28it/s]


Window 0: 100%|██████████| 1/1 [00:01<00:00,  1.37s/it]
2/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R022.tif, 6 files:  40%|████      | 2/5 [00:01<00:02,  1.37it/s]




Window 1: 100%|██████████| 2/2 [00:09<00:00,  4.78s/it]
3/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R008.tif, 6 files:  60%|██████    | 3/5 [00:11<00:09,  4.61s/it]




Window 0: 100%|██████████| 2/2 [00:09<00:00,  4.74s/it]
4/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R065.tif, 12 files:  80%|████████  | 4/5 [00:21<00:06,  6.61s/it]








Window 1: 100%|██████████| 4/4 [00:16<00:00,  4.22s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R108.tif, 13 files: 100%|██████████| 5/5 [00:39<00:00, 10.50s/it]








Window 1: 100%|██████████| 4/4 [00:17<00:00,  4.42s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R108.tif, 13 files: 100%|██████████| 5/5 [00:57<00:00, 11.59s/it]

/data/gh/tutorial/tiles/s2l2a_summer2019_B08/P50/
[2022-06-29 18:12:26] mosaic: /data/gh/tutorial/tiles/s2l2a_summer2019_B08/P50/
[2022-06-29 18:12:26] number of orbits: 5


1/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R051.tif, 1 files:  20%|██        | 1/5 [00:00<00:00, 12557.80it/s]


Window 0: 100%|██████████| 1/1 [00:01<00:00,  1.39s/it]
2/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R022.tif, 6 files:  40%|████      | 2/5 [00:01<00:02,  1.35it/s]




Window 1: 100%|██████████| 2/2 [00:09<00:00,  4.78s/it]
3/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R008.tif, 6 files:  60%|██████    | 3/5 [00:11<00:09,  4.63s/it]




Window 0: 100%|██████████| 2/2 [00:09<00:00,  4.69s/it]
4/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R065.tif, 12 files:  80%|████████  | 4/5 [00:21<00:06,  6.58s/it]








Window 3: 100%|██████████| 4/4 [00:17<00:00,  4.28s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R108.tif, 13 files: 100%|██████████| 5/5 [00:39<00:00, 10.57s/it]








Window 1: 100%|██████████| 4/4 [00:17<00:00,  4.45s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R108.tif, 13 files: 100%|██████████| 5/5 [00:58<00:00, 11.65s/it]

/data/gh/tutorial/tiles/s2l2a_summer2019_B08/P75/
[2022-06-29 18:13:24] mosaic: /data/gh/tutorial/tiles/s2l2a_summer2019_B08/P75/
[2022-06-29 18:13:24] number of orbits: 5


1/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R051.tif, 1 files:  20%|██        | 1/5 [00:00<00:00, 12192.74it/s]


Window 0: 100%|██████████| 1/1 [00:01<00:00,  1.35s/it]
2/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R022.tif, 6 files:  40%|████      | 2/5 [00:01<00:02,  1.38it/s]




Window 1: 100%|██████████| 2/2 [00:09<00:00,  4.81s/it]
3/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R008.tif, 6 files:  60%|██████    | 3/5 [00:11<00:09,  4.63s/it]




Window 0: 100%|██████████| 2/2 [00:09<00:00,  4.82s/it]
4/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R065.tif, 12 files:  80%|████████  | 4/5 [00:21<00:06,  6.68s/it]








Window 1: 100%|██████████| 4/4 [00:17<00:00,  4.29s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R108.tif, 13 files: 100%|██████████| 5/5 [00:39<00:00, 10.65s/it]








Window 1: 100%|██████████| 4/4 [00:17<00:00,  4.46s/it]
5/5 /data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R108.tif, 13 files: 100%|██████████| 5/5 [00:58<00:00, 11.72s/it]

9.1.2. 3. Final stitching¶

First we need to find the extent of the final mosaic.

[18]:

left, bottom = 4154000, 3145000
right, top = 4564000, 3446000
extent = [left, bottom, right, top]

Then we can run the stitching method for each statistic.

[19]:

for p in ['CNT','P25','P50','P75']:
    folder_input = f'/data/gh/tutorial/orbits/s2l2a_summer2019_B08/{p}/'
    file_output = f'/data/gh/tutorial/final/s2l2a_summer2019_B08_{p}.tif'
    print(folder_input)
    mosaic_maker.mosaic_final_tiled(folder_input, file_output, extent=extent)

/data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/
[2022-06-29 18:14:40] /data/gh/tutorial/final/s2l2a_summer2019_B08_CNT.tif
[2022-06-29 18:14:40] Downloading orbit files ...

/data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R051.tif: : 4it [00:00, 3013.69it/s]

[2022-06-29 18:14:40]

Done.
[2022-06-29 18:14:40]   mosaic


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/CNT/R051.tif: : 5it [00:00, 895.03it/s]




Window 2: 100%|██████████| 2/2 [00:09<00:00,  4.87s/it]

[2022-06-29 18:14:50] None
/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/
[2022-06-29 18:14:50] /data/gh/tutorial/final/s2l2a_summer2019_B08_P25.tif
[2022-06-29 18:14:50] Downloading orbit files ...


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R051.tif: : 4it [00:00, 3110.35it/s]

[2022-06-29 18:14:50]

Done.
[2022-06-29 18:14:50]   mosaic


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P25/R051.tif: : 5it [00:00, 769.03it/s]




Window 2: 100%|██████████| 2/2 [00:08<00:00,  4.46s/it]

[2022-06-29 18:14:59] None
/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/
[2022-06-29 18:14:59] /data/gh/tutorial/final/s2l2a_summer2019_B08_P50.tif
[2022-06-29 18:14:59] Downloading orbit files ...


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R051.tif: : 4it [00:00, 2367.66it/s]

[2022-06-29 18:14:59]

Done.
[2022-06-29 18:14:59]   mosaic


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P50/R051.tif: : 5it [00:00, 664.10it/s]



Window 1: 100%|██████████| 2/2 [00:08<00:00,  4.38s/it]

[2022-06-29 18:15:08] None
/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/
[2022-06-29 18:15:08] /data/gh/tutorial/final/s2l2a_summer2019_B08_P75.tif
[2022-06-29 18:15:08] Downloading orbit files ...


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R051.tif: : 4it [00:00, 2264.74it/s]

[2022-06-29 18:15:08]

Done.
[2022-06-29 18:15:08]   mosaic


/data/gh/tutorial/orbits/s2l2a_summer2019_B08/P75/R051.tif: : 5it [00:00, 769.15it/s]



Window 1: 100%|██████████| 2/2 [00:08<00:00,  4.47s/it]

[2022-06-29 18:15:17] None