venerdì 3 febbraio 2023

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




martedì 24 gennaio 2023

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

giovedì 19 gennaio 2023

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)




lunedì 16 gennaio 2023

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


mercoledì 11 gennaio 2023

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

mercoledì 14 dicembre 2022

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

martedì 6 dicembre 2022

Subset geografico con SNAP

Ci ho messo un po' a capire come funziona....

Per effettuare il subset via graph in SNAP si devono selezionare le bande che si intendono utilizzare e la finestra geografica. Per farlo non in modo interattivo si deve inserire un poligono in formato WKT e georiferito in WGS 84 (anche se le immagini sono in altro sistema di riferimento) ..partendo da uno shapefile si puo' utilizzare il servizio on line di conversione e poi fare update (se tutto e' andato bene compare un poligono giallo all'interno del perimetro rosso dell'immagine da tagliare




Debugger integrato ESP32S3

Aggiornamento In realta' il Jtag USB funziona anche sui moduli cinesi Il problema risiede  nell'ID USB della porta Jtag. Nel modulo...