Download and process Sentinel 5P data

 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()

Per effettuare una query geografica si puo' usare il motore di postgis indicando il punto geografico desiderato ed un buffer..con questa query i dati dello stesso giorno vengono mediati

 

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)