Query Sentinel-2 images¶
In this notebook we will show how to query Sentinel-2 images over an area between two given dates. We will use the sentinelsat package to obtain the Sentinel-2 product names and we will read those from the public GCP bucket (gs://gcp-public-data-sentinel-2/tiles).
Set the env variables to be able to read from the Google bucket. This will incur in reading costs.
Install package with Google dependecies¶
This is needed to read image from S2 bucket and to query to Copernicus SciHub
pip install georeader-spaceml fsspec gcsfs sentinelsat
import os
import datetime
from shapely.geometry import box
import geopandas as gpd
from georeader.readers import S2_SAFE_reader
# Donwload key from next line link to access the buckets and requester pays requests to public bucket (this is needed to query Sentinel-2 data)
# This is required to do advaced operations in the GCP bucket
# os.environ["GOOGLE_APPLICATION_CREDENTIALS"] = "path/to/file.json"
# os.environ["GS_USER_PROJECT"] = "project-name"
# S2_SAFE_reader.DEFAULT_REQUESTER_PAYS=True
os.environ["GS_NO_SIGN_REQUEST"] = "YES"
Create a geojson file using geojson.io
import geopandas as gpd
# aoi = gpd.read_file("/home/gonzalo/Downloads/sagunt_small.geojson")
aoi = gpd.read_file("/home/gonzalo/Downloads/liria.geojson")
# aoi.explore()
Step 1: Select dates and area of interest to read¶
from zoneinfo import ZoneInfo
tz = ZoneInfo("Europe/Madrid")
polygon_read = aoi.unary_union
crs_polygon = aoi.crs
datetime_str = "2019-09-28"
date_of_interest = datetime.datetime.strptime(datetime_str, "%Y-%m-%d").replace(tzinfo=tz)
date_start_search = date_of_interest - datetime.timedelta(days=10)
date_end_search = date_of_interest + datetime.timedelta(days=6)
print(f"Querying images between {date_start_search} and {date_end_search}\nArea: {polygon_read}")
Step 2: Query the products¶
from georeader.readers import scihubcopernicus_query
# from sentinelsat.sentinel import SentinelAPI
# # 'https://scihub.copernicus.eu/apihub'
# api = SentinelAPI('gonzmg88', "yyyyyyy", api_url='https://scihub.copernicus.eu/dhus/')
products_gpd = scihubcopernicus_query.query(polygon_read, date_start_search, date_end_search)
products_gpd[["overlappercentage","cloudcoverpercentage", "utcdatetime","localdatetime","solardatetime","solarday"]]
Show products queried¶
import folium
products_gpd["localdatetime_str"] = products_gpd["localdatetime"].dt.strftime("%Y-%m-%d %H:%M:%S")
m = products_gpd[["geometry","overlappercentage","cloudcoverpercentage", "localdatetime_str","solarday"]].explore("solarday",name="S2")
aoi.explore(m=m,name="AoI",color="red")
folium.LayerControl(collapsed=False).add_to(m)
m
Step 3: Read & plot the data¶
Here we will loop over the solar days querying and mosaicking the images over the AoI.
%%time
import numpy as np
from georeader import read
from georeader import mosaic
from georeader import window_utils
import rasterio.plot as rstplt
import matplotlib.pyplot as plt
for day, products_gpd_day in products_gpd.groupby("solarday"):
products_read = products_gpd_day.index
print(f"Selected {products_read.shape[0]} products for solar day {day}")
s2objs = []
for product in products_read:
s2_safe_folder = S2_SAFE_reader.s2_public_bucket_path(product+".SAFE", check_exists=False)
s2objs.append(S2_SAFE_reader.s2loader(s2_safe_folder, out_res=10, bands=["B04", "B03", "B02"]))
polygon_read_dst_crs = window_utils.polygon_to_crs(polygon_read, crs_polygon=crs_polygon, dst_crs=s2objs[0].crs)
data_memory = mosaic.spatial_mosaic(s2objs, polygon=polygon_read_dst_crs, dst_crs= s2objs[0].crs)
print(repr(data_memory))
fig, ax = plt.subplots(1,1,figsize=(8,8))
rstplt.show(np.clip((data_memory.values)/3_000,0,1), transform=data_memory.transform, ax=ax)
ax.set_title(day)
plt.show(fig)
plt.close(fig)
In the last day there's a missing image (for that reason we have a black square in the top left corner). We inspect the available products in the cell bellow.
import geopandas as gpd
pol_mosaic = data_memory.footprint(crs=crs_polygon)
print(f"Overlap {pol_mosaic.intersection(polygon_read).area / polygon_read.area*100:.2f}%")
footprint_downloaded = gpd.GeoDataFrame({"geometry": [pol_mosaic],
"title": ["footprint mosaic"],
"group": ["footprint mosaic"]},
crs=crs_polygon)
m = footprint_downloaded.explore(name="footprint",color="green")
m = aoi.explore(m=m,name="AoI",color="red")
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},
}