Overlapping S2 and Proba-V
Read overlapping images from Proba-V and Sentinel-2¶
Given a Sentinel-2 GeoTIFF file from the CloudSEN12 dataset and a Proba-V image we will show how to read the overlapping regions of those two using the georeader package.
The Sentinel-2 image from CloudSEN12 that we will use in this tutorial is: cloudsen12/high/point_0317/20190207T172509_20190207T173213_T13QGG/S2L1C.tif. In order to proceed you need to download the image and store it in the same directory as this tutorial.
Install package with Proba-V dependecies¶
pip install git+https://github.com/spaceml-org/georeader#egg=georeader[probav]
Download the Proba-V image
from georeader.readers import download_pv_product
link_pv_product = "https://www.vito-eodata.be/PDF/datapool/Free_Data/PROBA-V_100m/S1_TOA_100_m_C1/2019/2/9/PV_S1_TOA-20190209_100M_V101/PROBAV_S1_TOA_X07Y05_20190209_100M_V101.HDF5"
filename_down = download_pv_product.download_product(link_pv_product)
filename_down
Inspect Proba-V products¶
Prova-V images are very large (10,080x10,080 pixels). Therefore we should avoid read them all in memory if possible.
%%time
from georeader.readers import probav_image_operational
toa_reader = probav_image_operational.ProbaVRadiometry(filename_down, level_name="LEVEL3")
toa_reader
Warning!!¶
To read the Proba-V image content the h5py package requires an special compression. If installed from conda this compressor is not available. Therefore, re-install h5py from pip if the following cell fails!
pip install h5py --no-deps --ignore-installed
toa_reader.assert_can_be_read()
# If this raises an error reinstall h5py with pip and reset the notebook:
# pip install h5py --no-deps --ignore-installed
Read Sentinel-2 file of the cloudSEN12 dataset¶
from georeader.rasterio_reader import RasterioReader
s2reader = RasterioReader("S2L1C.tif")
s2reader
Band names are stored in the descriptions attribute
s2reader.descriptions
Set the reader to read the BLUE, RED NIR and SWIR bands of Sentinel-2 that are equivalent to the bands of Proba-V.
s2reader.set_indexes_by_name(["B2","B4", "B8","B11"])
s2reader
# Read image in memory
s2reader_memory = s2reader.load()
s2reader_memory.values
Read Proba-V image within the bounds of the Sentinel-2¶
%%time
from georeader import read
tile_read = read.read_from_bounds(toa_reader, bounds=s2reader.bounds, crs_bounds=s2reader.crs,
pad_add=(50, 50))
tile_read
tile_read_memory = tile_read.load()
tile_read_memory.values
%%time
import numpy as np
import matplotlib.pyplot as plt
# SWIR, NIR, RED composite
rgb_show = np.transpose(tile_read.values[:0:-1],(1,2,0))
plt.imshow(rgb_show)
plt.imshow(np.clip(s2reader_memory.values[:0:-1]/10_000,0,1).transpose((1,2,0)))
Reproject Proba-V and Sentinel-2¶
First we will reproject the Proba-V image to the UTM crs of Sentinel-2 at a 100m resolution
probav_100m_utm = read.read_reproject(toa_reader, bounds=s2reader.bounds, dst_crs=s2reader.crs,
resolution_dst_crs=100)
probav_100m_utm
Secondly we will reproject Sentinel-2 to 100m using an anti-aliasing filter to avoid artifacts.
s2reader_memory.set_dtype(np.float32) # Convert to float to avoid resampling errors!
s2_100m = read.resize(s2reader_memory, resolution_dst=100, anti_aliasing=True)
s2_100m
fig, ax = plt.subplots(2,3,figsize=(18,12))
ax[0,0].imshow(probav_100m_utm.values[:0:-1].transpose((1,2,0)))
ax[0,0].set_title("Proba-V 100m (UTM grid) SWIR/NIR/RED")
ax[0,1].imshow(np.clip(s2_100m.values[:0:-1].transpose((1,2,0))/10_000,0,1))
ax[0,1].set_title("Sentinel-2 100m SWIR/NIR/RED")
ax[0,2].imshow(np.clip(s2reader_memory.values[:0:-1].transpose((1,2,0))/10_000,0,1))
ax[0,2].set_title("Sentinel-2 10m SWIR/NIR/RED")
ax[1,0].imshow(probav_100m_utm.values[-2::-1].transpose((1,2,0)))
ax[1,0].set_title("Proba-V 100m (UTM grid) NIR/RED/BLUE")
ax[1,1].imshow(np.clip(s2_100m.values[-2::-1].transpose((1,2,0))/10_000,0,1))
ax[1,1].set_title("Sentinel-2 100m NIR/RED/BLUE")
ax[1,2].imshow(np.clip(s2reader_memory.values[-2::-1].transpose((1,2,0))/10_000,0,1))
ax[1,2].set_title("Sentinel-2 10m NIR/RED/BLUE")
Licence¶
The georeader package is published under a GNU Lesser GPL v3 licence
If you find this work useful please cite:
@article{portales-julia_global_2023,
title = {Global flood extent segmentation in optical satellite images},
volume = {13},
issn = {2045-2322},
doi = {10.1038/s41598-023-47595-7},
number = {1},
urldate = {2023-11-30},
journal = {Scientific Reports},
author = {Portalés-Julià, Enrique and Mateo-García, Gonzalo and Purcell, Cormac and Gómez-Chova, Luis},
month = nov,
year = {2023},
pages = {20316},
}