Frecce tricolori Firenze 11 febbraio 2023

 

zona Firenze Sud

Sentinel 5 NetCDF CH4

Ho provato a scaricare i dati del Sentinel 5, un satellite specializzato nel rilevamento di gas atmosferici (CH4, CO, formaldeide, O3,...)

I dati si possono scaricare da SciHub (https://s5phub.copernicus.eu/dhus/#/home) e sono distribuiti in formato NetCDF. Le dimensione del pixel sono di circa 4x7 Km

il dato e' da intendersi come in attesa di validazione (Sentinel 5 e' in fase pre-operativa). Come indicato anche nelle note di Google Earth Engine non solo basta scartare i pixel con il valore di qa <0.5 (e sono tanti in ogni immagine). Il valore di q5 si estrae direttamente dal file netcdf

Timelapse Luglio 2022 CH4 

I dati possono essere visualizzati con Planoply ma ho preferito provare ad utilizzare le librerie di Python 

from netCDF4 import Dataset

import numpy as np

from mpl_toolkits.basemap import Basemap

import sys

file = sys.argv[1]

data = Dataset(file, mode='r') # read the data

print(type(data)) # print the type of the data

print(data.groups['PRODUCT']) # print the variables in the data

print (data.groups['PRODUCT'].variables.keys())

ch4 = data.groups['PRODUCT'].variables['methane_mixing_ratio'][0,:,:]

lons = data.groups['PRODUCT'].variables['longitude'][:][0,:,:]

lats = data.groups['PRODUCT'].variables['latitude'][:][0,:,:]

print (lons.shape)

print (lats.shape)

print (ch4.shape)

ch4_units = data.groups['PRODUCT'].variables['methane_mixing_ratio'].units

import matplotlib.pyplot as plt

from matplotlib.colors import LogNorm

from mpl_toolkits.basemap import Basemap

#lon_0 = lons.mean()

#lat_0 = lats.mean()

lon_0 = 11

lat_0 = 43

m = Basemap(width=1000000,height=1500000,

resolution='l',projection='stere',\

lat_ts=2,lat_0=lat_0,lon_0=lon_0)

xi, yi = m(lons, lats)

# Plot Data

cs = m.pcolor(xi,yi,np.squeeze(ch4),norm=LogNorm(), cmap='jet')

# Add Grid Lines

m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)

m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries

m.drawcoastlines()

m.drawstates()

m.drawcountries()

# Add Colorbar

cbar = m.colorbar(cs, location='bottom', pad="10%")

cbar.set_label(ch4_units)

# Add Title

plt.title(file)

plt.get_current_fig_manager().full_screen_toggle() # toggle fullscreen mode

#plt.show()

plt.savefig(file+".png",dpi=300)

una elaborazione puo' essere fatta anche Google Earth Engine

=========================

var collection = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_CH4')

  .select('CH4_column_volume_mixing_ratio_dry_air')

  .filterDate('2019-01-01', '2019-12-31');

  

var band_viz = {

  min: 1000,

  max: 2000,

  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']

};

var median = collection.reduce(ee.Reducer.median());

print(median)

Map.addLayer(median, band_viz, 'S5P CH4');

Map.setCenter(11.295553516712024, 43.976202446636606,10);

PCA e GLCM in Sentinel 1

Utilizzando le immagini GRDH di Sentinel 1 (il dato GRDH si differenzia dal SLC perche' nel primo si ha solo l'ampiezza e l'intesita' del segnale nelle due polarizzazioni VV e VH mentre il SLC si ha sia la componente reale che quella complessa) si possono utilizzare i livelli di risposta (scale di grigi nelle immagini) per applicare degli algoritmi di PCA (componenti principali) e GLCM (Gray Level Co occurence matrix) per analisi delle texture

Dissimilarity map da GLCM

Composit di componenti principali

Due esempi di elaborazione in PCA e GLCM pre e post evento di inondazione

Partendo dal dato GRDH i passi da seguire sono (https://www.youtube.com/watch?v=xoQ4NikdOq0)

• Subset
• Apply Orbit File
• Radiometric -> Calibrate
• Speckle Filtering
• Range Doppler Terrain Correction
• Raster/Data Conversion/Convert to dB
• Raster/Image Analysis/ Principal components
• Raster/Image Analysis/Texture Analysis/Gray Level Co occurence Matrix 

Nel caso di PCA nelle bande troveremo le mappe per ogni componente principale con in aggiunta la mappa di errore e  response

Nel caso di GCLM troveremo per sigma0 VV e VH la mappa di contrasto, dissimilarita', omogeneita', ASM, energia, entropia, media, varianza,correlazione

Sentinel 2 Super resolution

 Mediante il plugin di SNAP Sen2Res e' possibile ottenere una migliore risoluzione spaziale (per esempio per passare da 20 m/px a 10 m/px)

Una volta installato nel menu Optical compare la voce Sentinel 2 Super Resolution

Il calcolo e' estremamente impegnativo dal punto di vista computazionale. A titolo di esempio per effettuare la super resolution su una sola banda su una porzione inferiore al 20% di tutta una immagine ha richiesto al mio portatile piu' di un'ora. Questo perche' l'algoritmo conserva le informazioni spetttrali nelle immagini e non e' come fare un pancromatico su Landsat (dove l'informazione di alta risoluzione e' inserita nel canale luminosita' di un composit RGB)

Vertical displacement from phase in Sentinel 1

Un approccio alternativo alla creazione di un DEM da interferogramma e' quello di usare Phase to Displacement al posto di Phase to Elevation  (in Radar/Interferometric/Products)

Ho creato un inteferogramma basandosi su una immagine Sentinel 1 prima e dopo un evento basandosi su questo tutorial

https://www.youtube.com/watch?v=9__baWNmoJQ

Dem da Sentinel 1

 Per creare un DEM da dati Sentinel 1 si devono per prima cosa selezionare due immagini SLC che siano il piu' possibili coerenti (quindi con date di acquisizione ravvicinate e la stessa traiettoria ascendente o discendente),con una baseline ampia (per selezionare le immagini si puo' utilizzare Vertex di ASF Alaska) ed immagini che appartengono alla stessa path

In generale si devono anche evitare immagini che presentino zone di mare (o laghi) e per il possibile di deve evitare periodi dell'anno con presenza di umidita'/precipitazioni. Anche le zone con spiccate modifiche della vegetazione possono creare problemi a causa della non penetrazione della banda C del radar sul canopy

Le immagini (i file sono gli SLC IW) sono molto grandi e quindi il primo passo e' effettuare uno split 

Dall'interno di SNAP si puo' selezionare sono la subswath e la polarizzazione ma si puo' selezionare un sottoinsieme ancora piu' ridotto selezionando i bursts creando un file come il successivo e caricandolo in S1- TOPS Split. Se non si vuole usare il file si puo' fare anche dall'interfaccia utente cambiando la slide bar orizzontale al di sotto della scritta Bursts

<parameters>

    <subswath>IW2</subswath>

    <selectedPolarisations>VV</selectedPolarisations>

    <firstBurstIndex>9</firstBurstIndex>

    <lastBurstIndex>10</lastBurstIndex>

</parameters>

In seguito si applica la correzione di orbita su ogni split (Radar/Apply Orbit File)

Si crea quindi lo stack delle due acquisizioni con S1 Back Geocoding (togliere Mask out area with no elevation ed flaggare Output Deramp and Demod Phase

• Creiamo il file con nome Stack . Dopo questa fase in modo opzionale si puo' effettuare S1 Enhanced Spectral Diversity
• Usando il comandoRadar/Interferometric/Products/Inteferogram formation si ottiene come risultato file Stack_ifg
• Si procede al debusrst conS1 TOPS Deburst file Stack_ifg_deb
• Il passo successivo e 'Radar > Interferometric > Products/Topographic Phase Removal option con file risultato  Stack_ifg_deb_dinsar
• Si applica Multilooking cambiando a 6 il parametro Number of Range Looks file Stack_ifg_deb_dinsar_ML

• Fase di filtraggio Radar > Interferometric > Filtering/Goldstein Phase Filtering file Stack_ifg_deb_dinsar_ML_flt
• Export to SNAPSU TOPO Number tile 20 20
• Si entra nella directory in cui e' stato fatto l'export e si seleziona il file di Phase per esempio in polarizzazione VV Phase_ifg_VV_25Nov2022_31Dec2022.snaphu.hdr e si fa Radar/Interferometric/Unwrapping/snaphu-unwrapping
• in seguito si fa l'import selezionando come Source Product il file  Stack_ifg_deb_dinsar_ML_flt e come Unwrapped il file unwrapped nel folder snaphu. Si salva come unwrapped
• Radar/Inteferometric/Products/Phase to elevation scegliendo il parametro SRTM 1 Sec HGT risultato file unwrapped_dem
• Geometric > Terrain Correction > Range Doppler Terrain Correction. si sceglie il parametro SRTM 1 Sec HGT si puo' anche scegliere la proiezione indicando WGS84/UTM e tra selected sources si flagga DEM file unwrapped_dem_TC. Si puo' quindi esportare la banda come Geotiff

Layer stack con RasterIO

La libreria RasterIO permette di fare layer stacking ovvero creare dei composit 

In questo caso sono stati presi 3 indici NDVI della stessa area in diversi periodi dell'anno 

Il primo file nell'array corrispondera' al canale rosso ed i successibi verde e blu rispettivamente (tutti i layer devono avere le stesse dimensioni)

import rasterio

file_list = ['20220804_ndvi.tiff','20220913_ndvi.tiff','20221028_ndvi.tiff']

# Read metadata of first file

with rasterio.open(file_list[0]) as src0:

meta = src0.meta

# Update meta to reflect the number of layers

meta.update(count = len(file_list))

# Read each layer and write it to stack

with rasterio.open('stack.tif', 'w', **meta) as dst:

for id, layer in enumerate(file_list, start=1):

with rasterio.open(layer) as src1:

dst.write_band(id, src1.read(1))