venerdì 17 febbraio 2023

Eliminare dati a bassa coerenza da inteferogrammi Sentinel 1

Nel caso i dati di coerenza di un inteferogramma sentinel 1 siano bassi (tali da rendere inutili ulteriori elaborazioni) si puo' tentare di filtrare i dati di fase con bassa coerenza


 

Per fare cio' si seleziona la banda della fase e con tasto destro si selezione Properties. In seguito si edita Pixel Values Expression con una formula tipo quella sottostante

if coh_IW1_VV_12Dec2022_24Dec2022 > 0.4 then atan2(atan2(q_ifg_IW1_VV_12Dec2022_24Dec2022,i_ifg_IW1_VV_12Dec2022_24Dec2022))  else 0

in pratica tutti i pixel con coerenza inferiore a 0.4 vengono messi a zero
Si salva il prodotto e si puo' procedere con i passi successivi

lunedì 13 febbraio 2023

DEM e coerenza su Sentinel 1

 Dopo un po' di insuccessi nella creazione di DEM da dati Sentinel 1 ho ricercato in bibliografia le possibili cause

Per prima cosa dopo aver fatto l'interferogramma si deve selezionare l'immagine di coerenza e creare l'istogramma dal menu Analysis/Histogram (premendo refresh),

In generale per avere un DEM significativo il valore di coerenza di almeno il 70% dei punti deve essere almeno 0.4


Scarsi valori di coerenza sono imputabili a cause atmosferiche

Viene indicato da vari autori un offset tra 20 e 40 m nei dati di elevazione tra il dato reale ed il dato radar

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)




Red Edge con Sentinel 2

Volevo migliorare un po' quanto provato qui piu' che altro per avere una migliore risoluzione spaziale. Ho provato con Sentinel 2 (...