Per scaricare in automatico i dati Sentinel ho utilizzato il progeto CDSETool
ho creato un file shape in coordinate WGS86 4236 nel folder dello script
from cdsetool.query import query_features, shape_to_wkt
from cdsetool.credentials import Credentials
from cdsetool.download import download_features
from cdsetool.monitor import StatusMonitor
from datetime import date
#chiavi di ricerca per Sentinel 5
# dict_keys(['maxRecords', 'index', 'page', 'identifier', 'geometry', 'box', 'lon', 'lat', 'radius', 'startDate', 'completionDate',
# 'productIdentifier', 'productType', 'processingLevel', 'platform', 'instrument', 'orbitNumber', 'sensorMode',
# 'updated', 'publishedAfter', 'publishedBefore', 'sortParam', 'sortOrder', 'status', 'exactCount', 'processingBaseline'])
PASSWORD_SENTINEL_ESA = "xxxxxxxx"
USER_SENTINEL_ESA = "l.innocenti@xxxxxx"
PATH_WHERE_TO_SAVE = "/home/luca/netcdf/"
geometry = shape_to_wkt("toscana.shp")
print(geometry)
features = query_features(
"Sentinel5P",
{
"startDate": "2024-05-30",
"completionDate": date(2024, 5, 31),
"maxRecords": 5,
#"processingLevel": "S5PL2",
"productType": "L2__O3____",
# title_ids : ["32TPP"], solo per Sentinel 2
"geometry": geometry
},
)
list(
download_features(
features,
PATH_WHERE_TO_SAVE,
{
"concurrency": 4,
"monitor": StatusMonitor(),
"credentials": Credentials(USER_SENTINEL_ESA, PASSWORD_SENTINEL_ESA),
},
)
)
ATTENZIONE: i files scaricati nonostante riportino una estensione nc in realta' sono zippati. Si devono quindi prima decomprimere e poi trattare come NetCDF.nc)
Per visualizzare in modo speditivo i dati si puo' usare Panoply. Alrimenti si puo' usare anche SNAP (per vedere i dati georiferiti https://www.youtube.com/watch?v=G8tVNbdu8-A)
il file nc puo' essere letto tramite xarray
in sintesi ci sono diversi array, uno contiene le misure, uno la longitudine del pixel ed uno la latitudine del pixel
La cosa piu' comoda per l'analisi e' salvare i dati in un database postgis. Attenzione che devono essere gestiti i valori nan che possono essere presenti
#!/usr/bin/env python
# coding: utf-8
import xarray as xr
import math
import psycopg
from shapely.geometry import LineString
from shapely import wkb
conn = psycopg.connect("dbname=postgis_db user=postgres password=Tvlgal55 host=localhost")
curs = conn.cursor()
#CREATE TABLE ozono(geom geometry, valore FLOAT not NULL, data TIMESTAMP without time zone not NULL)'
#data5p = xr.open_dataset('./S5P_NRTI_L2__O3_____20240529T122201_20240529T122701_34335_03_020601_20240529T130706/S5P_NRTI_L2__O3_____20240529T122201_20240529T122701_34335_03_020601_20240529T130706.nc', group='PRODUCT')
data5p = xr.open_dataset('./S5P_NRTI_L2__O3_____20240530T120201_20240530T120701_34349_03_020601_20240530T124545/S5P_NRTI_L2__O3_____20240530T120201_20240530T120701_34349_03_020601_20240530T124545.nc', group='PRODUCT')
#print(data5p)
#print(data5p.dims)
dimx= data5p.dims['scanline']
dimy=data5p.dims['ground_pixel']
tempo = data5p['ozone_total_vertical_column']['time']
data= "2024-05-30 12:02:01"
#print(dimx)
#print(dimy)
#print(tempo)
for x in range(0,dimx):
for y in range(0,dimy):
lon = str(data5p.longitude.values[0,x,y])
lat = str(data5p.latitude.values[0,x,y])
valore = data5p.ozone_total_vertical_column.values[0,x,y]
#print(str(data5p.longitude.values[0, x,y])+";"+str(data5p.latitude.values[0, x,y])+";"+str(data5p.ozone_total_vertical_column.values[0, x,y]))
if not(math.isnan(valore)):
curs.execute("INSERT INTO ozono (geom, valore,data) VALUES (ST_SetSRID(ST_MakePoint("+lon+","+lat+"), 4326),"+str(valore)+",'"+data+"')");
conn.commit()
SELECT avg(valore) as median FROM ozono WHERE ST_DistanceSphere(geom, ST_MakePoint(11.1,43.1)) <= 10000 GROUP BY data
Si puo' usare QGis per vestire i dati con un stile ed esportando il file .SLD
La vestizione puo' essere importata in Geoserver
Il prodotto finale puo' essere distribuito con Openlayers su web
per indicizzare un campo spaziale in postgis si puo' usare la sintassi
CREATE INDEX idx_ozono ON ozono USING gist (geom);
i dati in geoserver possono essere filtrati tramite il campo CQL basandosi su un box oppure tramite un buffer attorno ad un punto
BBOX(geom, 10, 42.5, 12, 44.5)
DWithin(geom,POINT(11.38 43.3),140000,meters)
