Visualizzazione post con etichetta Telerilevamento. Mostra tutti i post
Visualizzazione post con etichetta Telerilevamento. Mostra tutti i post

venerdì 21 aprile 2023

Earth Engine Sentinel 5P time series NO2

 Due script per scaricare serie tempo dei dati Sentinel 5P da Google Earth Engine dato un punto geografico




Tutti i dati

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

var point = ee.Geometry.Point([11.097849, 43.793234]);
var sentinel = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2');
var sentinelLST = sentinel.filterBounds(point)
.filterDate('2019-01-01', '2022-12-31')
.select('NO2_column_number_density');

sentinelLST = sentinelLST.map(function(img){
var date = img.get('system:time_start');
return img.multiply(100000).set('system_time_start', date);
});

var createTS = function(img){
var date = img.get('system_time_start');
var value = img.reduceRegion(ee.Reducer.mean(), point).get('NO2_column_number_density');
var ft = ee.Feature(null, {'system:time_start': date,
'date': ee.Date(date).format('Y/M/d'),
'value': value});
return ft;
};

var TS = sentinelLST.map(createTS);

var graph = ui.Chart.feature.byFeature(TS, 'system:time_start', 'value');

print(graph.setChartType("Sentinel 5P")
.setOptions({vAxis: {title: 'NO2'},
hAxis: {title: 'Date'}}));


Export.table.toDrive({collection: TS, selectors: 'date, value'});

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

Medie Giornaliere

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

var collection = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')
  .select('tropospheric_NO2_column_number_density')
var daily_data = ee.ImageCollection(ee.List.sequence(2019,2019).map(function(year){
  var date1 = ee.Date.fromYMD(year,1,1)
  var date2 = date1.advance(1,'year')
  //Calculate the number of days per year
  var doy = date2.difference(date1,'day')
  var doyList = ee.List.sequence(1,doy)
  //Daily image mean synthesis using doy
  var day_imgs = doyList.map(function(doy){
    doy = ee.Number(doy)
    var temp_date = date1.advance(doy.subtract(1),"day")
    var temp_img = collection.filterDate(temp_date,temp_date.advance(1,'day'))
    return temp_img.mean().set("system:time_start",temp_date.millis())
  })
  return day_imgs
}).flatten())
Map.addLayer(daily_data)
var chart = ui.Chart.image.series({
imageCollection:daily_data,
region:roi,
reducer:ee.Reducer.mean(),
scale:1113.2,
// xProperty:,
})
print(chart)

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


Medie Mensili

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


var point = ee.Geometry.Point([11.097849, 43.793234]);
var dataset = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
          .filterBounds(point)
          .filterDate('2019-01-01', '2022-12-31')
          .select('NO2_column_number_density');
var months = ee.List.sequence(1, 12);
var start_year = 2019;
var start_date = '2019-01-01';
var end_year = 2022;
var end_date = '2022-12-31';
var years = ee.List.sequence( start_year, end_year);

var byMonthYear =  ee.FeatureCollection(
  years.map(function (y) {
    return months.map(function(m){
          var w = dataset.filter(ee.Filter.calendarRange(y, y, 'year'))
                    .filter(ee.Filter.calendarRange(m, m, 'month'))
                    .mean();
      var pointMean = w.reduceRegion({reducer:ee.Reducer.first(), geometry:point,scale:1000});  
      return ee.Feature(null).set("value",pointMean.get("NO2_column_number_density")).set("year",y).set("month",m);
    })
  }).flatten()
);

print("feature collection",byMonthYear);

Export.table.toDrive({collection:byMonthYear,description:"O3_signa"})

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


mercoledì 8 marzo 2023

Formaldeide e metano da Sentinel 5 con Google Earth Engine

 Usare Google Earth Engine per Sentinel 5 e' comodo perche' i pixel sono stati ricalcolati in modo omogeneo (nei file netcdf originali i pixel hanno dimensioni differenti a secondo dell'orbita)

I valori di CH4 sono espressi in ppm ed il range di valore e' estremamente stretto (da 1860 a 1890 ppm nell'immagine ) con la dimensione del pixel di 1113 m. I dati sono disponibili sono per le aree di fondovalle




I dati di formaldeide sono espressi in mol/m2 ed il range di valori varia tra 0.00007 e 0.00014 con pixel di 1113 m. 




var palettes = require('users/gena/packages:palettes');
  
var palette_r = palettes.colorbrewer.Accent[5];
var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_HCHO')
  .select('tropospheric_HCHO_column_number_density')
  .filterDate('2019-01-01', '2019-12-31');
var median = collection.reduce(ee.Reducer.mean());
var band_viz = {
  min: 0.00007,
  max: 0.00014,
  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};
Map.addLayer(collection.mean(), band_viz, 'S5P HCHO');
Map.setCenter(11.295553516712024, 43.976202446636606,9);
Map.style().set('cursor', 'crosshair');
var header = ui.Label('HCHO', {fontWeight: 'bold', fontSize: '18px'});
var toolPanel = ui.Panel([header], 'flow', {width: '400px'});
Map.onClick(function(coords) {
  // Create or update the location label (the second widget in the panel)
  var location = 'lon: ' + coords.lon.toFixed(4) + ' ' +
                 'lat: ' + coords.lat.toFixed(4);
  var click_point = ee.Geometry.Point(coords.lon, coords.lat);
  var CH4Value = median.reduceRegion(ee.Reducer.first(), click_point, 90).evaluate(function(val){
    console.log(val);
    var CH4Text = 'CH4: ' + val;
    toolPanel.widgets().set(2, ui.Label(CH4Text));
  });
  
  toolPanel.widgets().set(1, ui.Label(location));
  // Add a red dot to the map where the user clicked.
  Map.layers().set(1, ui.Map.Layer(click_point, {color: 'FF0000'}));
});
// Add the panel to the ui.root.
ui.root.add(toolPanel);

giovedì 23 febbraio 2023

Composit Sentinel 2

Per creare un composit da Snap con dati Sentinel 2 si deve fare prima un subset (Raster/Subset) delle tre bande richieste. Se la risoluzione delle bande e' differente si deve successivamente rendere la grandezza dei pixel omogenea con Geometric/Resampling

    431

A questo punto Open RGB Image Window....in Colour Management si puo' impostare lo stretch


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

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

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




giovedì 24 novembre 2022

Grafico inteferogramma Snap

Si puo' automatizzare il lavoro visto nel precedente post medianto graph builder di SNAP

Come indicato da questo post (https://creodias.eu/-/interferometry-processing-example-) non e' possibile stabilire una pipeline ma si deve dividere il lavoro 3 step successivi

Step 1
(lo step 1 deve essere effettuato per ognuna delle due immagini di input)
Step 2

Step 3





mercoledì 23 novembre 2022

Interferogramma Sentinel 1

Ho provata ad applicare il tutorial per creare un interoferogramma da dati radar Sentinel 1 disponibile 

usando questi due acquisizioni del 10 marzo 3 17 novembre 2022. I dati sono stati selezionati tramite l'applicativo visto nel precedente post   

S1A_IW_SLC__1SDV_20220310T171440_20220310T171508_042262_05097D_2870.zip
S1A_IW_SLC__1SDV_20221117T171451_20221117T171518_045937_057F23_B553.zip



In tutto sono stati eseguiti 7 step 

  • Coregistration with ESD
  • Inteferogram formation
  • Deburst
  • Topographic Phase
  • Multilooking
  • Filtering 
  • Geocoding

Ho usato un normale laptop con un I5 6 gen ma con ampio spazio disco disponibile e 16Gb di ram

Si osservano delle chiare frange alle coordinate lat 43.90 lon 11.05 in localita' di Bagnolo di Prato

Utilizzando l'applicativo https://egms.land.copernicus.eu/ si puo' verificare che si tratta di un movimento reale e non di una mia elaborazione errata




Ricerca baseline per DEM Sentinel 1

Per cercare coppie di immagini radar Sentinel 1 per elaborare interferogrammi un sistema comodo e' utilizzare https://search.asf.alaska.edu/#/


Per la ricerca il metodo piu' comodo e' selezionare

Search Type : Geographic Search
Dataset : Sentinel 1
Draw Polygon

Dopo il Search si clicca su Baseline e si ha la selezione delle immagini.


Successivamente si possono scaricare i dati scaricando il file Python 


 

domenica 20 novembre 2022

Snap GPT Resample, Subset, MathBand

Resample a 10m, Subset su finestra geografica, calcolo NDVI


 


<graph id="Graph">
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <useAdvancedOptions>false</useAdvancedOptions>
      <file>${ingresso}</file>
      <copyMetadata>true</copyMetadata>
      <bandNames/>
      <pixelRegion>0,0,10980,10980</pixelRegion>
      <maskNames/>
    </parameters>
  </node>
  <node id="Resample">
    <operator>Resample</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <referenceBand/>
      <targetWidth>10980</targetWidth>
      <targetHeight>10980</targetHeight>
      <targetResolution/>
      <upsampling>Nearest</upsampling>
      <downsampling>First</downsampling>
      <flagDownsampling>First</flagDownsampling>
      <resamplingPreset/>
      <bandResamplings/>
      <resampleOnPyramidLevels>true</resampleOnPyramidLevels>
    </parameters>
  </node>
  <node id="Subset">
    <operator>Subset</operator>
    <sources>
      <sourceProduct refid="Resample"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <sourceBands>B4,B8</sourceBands>
      <tiePointGrids/>
      <region>0,0,0,0</region>
      <referenceBand/>
      <geoRegion>POLYGON ((11.141773223876953 43.80547332763672, 11.368977546691895 43.80547332763672, 11.368977546691895 43.664669036865234, 11.141773223876953 43.664669036865234, 11.141773223876953 43.80547332763672, 11.141773223876953 43.80547332763672))</geoRegion>
      <subSamplingX>1</subSamplingX>
      <subSamplingY>1</subSamplingY>
      <fullSwath>false</fullSwath>
      <copyMetadata>true</copyMetadata>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Subset"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
        <targetBand>
          <name>newBand</name>
          <type>float32</type>
          <expression>(B8 - B4) / (B8 + B4)</expression>
          <description/>
          <unit/>
          <noDataValue>0.0</noDataValue>
        </targetBand>
      </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <file>/home/luca/sentinel/Subset_resampled_BandMath.tif</file>
      <formatName>GeoTIFF</formatName>
    </parameters>
  </node>
  <applicationData id="Presentation">
    <Description/>
    <node id="Read">
            <displayPosition x="37.0" y="134.0"/>
    </node>
    <node id="Resample">
      <displayPosition x="148.0" y="141.0"/>
    </node>
    <node id="Subset">
      <displayPosition x="264.0" y="142.0"/>
    </node>
    <node id="BandMaths">
      <displayPosition x="359.0" y="144.0"/>
    </node>
    <node id="Write">
            <displayPosition x="455.0" y="135.0"/>
    </node>
  </applicationData>
</graph>


venerdì 18 novembre 2022

Landsat Google Cloud Bucket

Oltre alle immagini Sentinel e' possibile scaricare da Google  anche i dati Landsat Collection 1



Si trovano su un bucket pubblico

 https://cloud.google.com/storage/docs/public-datasets/landsat

il contenuto del folder indicizzato dal file index.csv.gz

il contenuto dell'indice (dimensione di oltre 1 Gb)contiene i seguenti campi

SCENE_ID,
PRODUCT_ID,
SPACECRAFT_ID,
SENSOR_ID,
DATE_ACQUIRED,
COLLECTION_NUMBER,
COLLECTION_CATEGORY,
SENSING_TIME,
DATA_TYPE,
WRS_PATH,
WRS_ROW,
CLOUD_COVER,
NORTH_LAT,
SOUTH_LAT,
WEST_LON,
EAST_LON,
TOTAL_SIZE,
BASE_URL

non e' disponibile un sistema di interrogazione non usando BigQuery...in pratica si scariva il csv e lo si elabora. Un metodo sbrigativo puo' essere l'utilizzo di csvq , lento ma che permette una sintassi in stile SQL. Per esempio

csvq "select * from index where WRS_PATH=192 and WRS_ROW=30

tramite la base URL si puo' scaricare il prodotto come da esempio successivo utilizzando le gsutil    

gsutil -m cp -r gs://gcp-public-data-landsat/LC08/01/044/034/LC08_L1GT_044034_20130330_20170310_01_T2/ .


Elaborazione batch con SNAP GPT headless

 E' possibile utilizzare SNAP in modalita' headless installando da remoto su un server il normale pacchetto di SNAP. A seconda del fatto che sia presente o meno una sessione di interfaccia grafica o meno viene eseguito un tipo diverso di installer

Per eseguire una catena di comandi in SNAP si puo' in maniera interattiva creare un grafico con Graph Builder e salvare il grafico in formato xml



questo e' un esempio di grafico salvato da SNAP per il calcolo NDVI tramite BandMath

=============================================================
<graph id="Graph">
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <useAdvancedOptions>false</useAdvancedOptions>
      <file>/home/luca/transi/S2B_MSIL2A_20220811T100559_N0400_R022_T32TPP_20220811T162101.zip</file>
      <copyMetadata>true</copyMetadata>
      <bandNames/>
      <pixelRegion>0,0,10980,10980</pixelRegion>
      <maskNames/>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
        <targetBand>
          <name>NDVI</name>
          <type>float32</type>
          <expression>(B8 - B4) / (B8 + B4)</expression>
          <description/>
          <unit/>
          <noDataValue>0.0</noDataValue>
        </targetBand>
      </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <file>/home/luca/transi/BandMath.tif</file>
      <formatName>GeoTIFF</formatName>
    </parameters>
  </node>
  <applicationData id="Presentation">
    <Description/>
    <node id="Read">
            <displayPosition x="37.0" y="134.0"/>
    </node>
    <node id="BandMaths">
      <displayPosition x="184.0" y="136.0"/>
    </node>
    <node id="Write">
            <displayPosition x="455.0" y="135.0"/>
    </node>
  </applicationData>
</graph>

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

per renderlo utilizzabile in modalita' batch si devono modificarlo per inserire delle variabili
Le variabili sono nel formato ${nome_variabile}
Nel codice sottostante e' stato reso variabile il file di input
=============================================================
<graph id="Graph">
  <version>1.0</version>
  <node id="Read">
    <operator>Read</operator>
    <sources/>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <useAdvancedOptions>false</useAdvancedOptions>
      <file>${ingresso}</file>
      <copyMetadata>true</copyMetadata>
      <bandNames/>
      <pixelRegion>0,0,10980,10980</pixelRegion>
      <maskNames/>
    </parameters>
  </node>
  <node id="BandMaths">
    <operator>BandMaths</operator>
    <sources>
      <sourceProduct refid="Read"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <targetBands>
        <targetBand>
          <name>NDVI</name>
          <type>float32</type>
          <expression>(B8 - B4) / (B8 + B4)</expression>
          <description/>
          <unit/>
          <noDataValue>0.0</noDataValue>
        </targetBand>
      </targetBands>
      <variables/>
    </parameters>
  </node>
  <node id="Write">
    <operator>Write</operator>
    <sources>
      <sourceProduct refid="BandMaths"/>
    </sources>
    <parameters class="com.bc.ceres.binding.dom.XppDomElement">
      <file>/home/luca/transi/ndvi.tif</file>
      <formatName>GeoTIFF</formatName>
    </parameters>
  </node>
  <applicationData id="Presentation">
    <Description/>
    <node id="Read">
            <displayPosition x="37.0" y="134.0"/>
    </node>
    <node id="BandMaths">
      <displayPosition x="184.0" y="136.0"/>
    </node>
    <node id="Write">
            <displayPosition x="455.0" y="135.0"/>
    </node>
  </applicationData>
</graph>
=============================================================

da linea di comando si puo' usare lo switch -P per popolare la variabile

-Pnome_variabile=valore_variabile

per lanciare il processo per esempio si puo' usare

/home/luca/snap/bin/gpt ndvi2.xml -Pingresso=/home/luca/transi/S2B_MSIL2A_20220811T100559_N0400_R022_T32TPP_20220811T162101.zip

giovedì 17 novembre 2022

Download script per Sentinel

Alcune volte ho avuto problemi nello scaricare dataset da SciHub a causa di connessioni interrotte a meta'

Per ovviare al problema mi sono scritto questo script bash che verifica la correttezza del file zip

#!/bin/bash
# yes | mv ./immagini/*.zip ./backup
#./dhusget_0.3.8.sh -u username -p password -m Sentinel-2 -t96 -o product -c 10.6,43.5:11.3,44 -n 1 -w 5 -O ./immagini

for file in `ls ./immagini/*.zip`
do
#echo $file
if unzip -t $file > /dev/null 2>&1
then
echo "Archivio $file OK"
#unzip $file
else
echo "Archivio $file corrotto"
rm -f $file
fi
done

Snap e Snappy

Si puo' automatizzare l'analisi delle immagini Sentinel, generalmente eseguite in SNAP, mediante il modulo python Snappy




Il problema e' che su Ubuntu 22.04 la procedura di compilazione automatica di Snappy fallisce ed anche in modo manuale non ho avuto successo

La via piu' breve e' quella di utilizzare un docker container. Quello aggiornato a Snap 8 non ha tutte le librerie Python gia' installato e quindi ho modificato il docker originale mediante il seguente Dockerfile

FROM mundialis/esa-snap:ubuntu
RUN apt-get update && apt-get -y install mc nano
RUN pip3 install --no-cache-dir numpy scipy pandas matplotlib

una volta scritto il Dockerfile si puo' effettuare la build del nuovo container con

docker build -t esa-snappy-numpy:0.1 .

per scambiare i file tra il container e il sistema operativo ho montato la home locale sulla root folder del container

 a questo punto si puo' procedere con una elaborazione per esempio un calcolo di NDVI. Si salva il seguente codice nel file ndvi.py in ./transi/ndvi.py 

Attenzione : nello stesso container e' possibile utilizzare anche il comando gpt per effettuare analisi batch usando il file di Graph Builder (creato tramite GUI all'esterno del container) /usr/local/snap/bin/gpt

il file e' stato preso da questo indirizzo GitHub. Un tutorial sull' uso di Snappy sono  qui e qui. Un esempio di preprocessing per Sentinel 1 e' disponibile qui

import sys
import numpy
from snappy import Product, ProductData, ProductIO, ProductUtils
if len(sys.argv) < 2:
    print ('Product file requires')
    sys.exit(1)
# input product & dimensions
input_product = ProductIO.readProduct(sys.argv[1])
width = input_product.getSceneRasterWidth()
height = input_product.getSceneRasterHeight()
product_name = input_product.getName()
# input product red & nir bands
red_band = input_product.getBand('B4')
nir_band = input_product.getBand('B8')
# output product (ndvi) & new band
output_product = Product('NDVI', 'NDVI', width, height)
ProductUtils.copyGeoCoding(input_product, output_product)
output_band = output_product.addBand('ndvi', ProductData.TYPE_FLOAT32)
# output writer
output_product_writer = ProductIO.getProductWriter('BEAM-DIMAP')
output_product.setProductWriter(output_product_writer)
output_product.writeHeader(product_name + '.ndvi.dim')
# compute & save ndvi line by line
red_row = numpy.zeros(width, dtype=numpy.float32)
nir_row = numpy.zeros(width, dtype=numpy.float32)
for y in range(height):
    red_row = red_band.readPixels(0, y, width, 1, red_row)
    nir_row = nir_band.readPixels(0, y, width, 1, nir_row)
    ndvi = (nir_row - red_row) / (nir_row + red_row)
    output_band.writePixels(0, y, width, 1, ndvi)
output_product.closeIO()

il comando readProduct legge lo zip completo scaricato da SciHub. Non si puo' modificare il nome del file

sudo docker run -it  -v /home/luca:/root esa-snappy-numpy:0.1 -c "python3.6 /root/transi/ndvi.py /root/transi/S2B_MSIL2A_20220811T100559_N0400_R022_T32TPP_20220811T162101.zip"

in uscita si ha come risultato un file leggibile in SNAP con l'elaborazione NDVI di tutto il granulo

il container si puo' usare anche in modo interattivo tramite bash digitando semplicemente
sudo docker run -it  -v /home/luca:/root esa-snappy-numpy:0.1

mercoledì 16 novembre 2022

Soil Moisture Sentinel 1 Google Earth Engine

 Visto che le immagini radar presenti su Google EarthEngine sono gia' radiometricamente e geometricamente corrette ho voluto provare a replicare su questa piattaforma quanto visto nel precedente post



la demo funzionante e' all'indirizzo sotto riportato

https://lucainnoc.users.earthengine.app/view/soil-moisture-sentinel-1

Questa e' la spvrapposizione tra la geologia e l'algoritmo dell'area del precedente post (Casole, Vicchio di Mugello). 

La mappa e' la differenza di umidita' tra il periodo primaverile ed estivo dell'anno 2021 con il colore blu che indica massima variazione negativa e il colore rosso massima variazione positiva


function maskS2clouds(image) {
  var qa = image.select('QA60');
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}

var imgVV = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .select('VV')
        .map(function(image) {
          var edge = image.lt(-30.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage);
        });
//seleziona solo ASCENDENTE
var asc = imgVV.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));

//seleziona il periodo primaverile ed estivo
var spring = ee.Filter.date('2021-03-01', '2021-04-20');
var summer = ee.Filter.date('2021-06-11', '2021-08-31');
//seleziona il periodo dal gennaio 2019 al dicembre 2021
//per ottenere i valori massimi e minimi
var full = ee.Filter.date('2019-01-01', '2021-12-31');
var sar_dataset_min = ee.Image.cat(asc.filter(full).min());
var sar_dataset_max = ee.Image.cat(asc.filter(full).max());
// calcola la variazione massima nel periodo di interesse 2019-2021
var sensitivity = sar_dataset_max.subtract(sar_dataset_min);

var primavera = ee.Image.cat(asc.filter(spring).mean());
var estate = ee.Image.cat(asc.filter(summer).mean());
var Mv_primavera = (primavera.subtract(sar_dataset_min)).divide(sensitivity);
var Mv_estate = (estate.subtract(sar_dataset_min)).divide(sensitivity);
var delta = Mv_estate.subtract(Mv_primavera);
//var Mv_estate = (estate - sar_dataset_min)/(sensitivity);
//print(Mv_primavera);
//estre i truecolor ottici per riferimento
var rgb_spring = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                  .filterDate('2018-03-01', '2018-04-20')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds);
var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};
var rgb_summer = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                  .filterDate('2018-06-11', '2018-08-31')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds);
var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

Map.setCenter(11.489351719442514,43.95792446517467, 14);
//Map.addLayer(primavera, {min: -25, max: 5}, 'Spring', true);
//Map.addLayer(estate, {min: -25, max: 5}, 'Summer', true);
Map.addLayer(rgb_spring.mean(), visualization, 'RGB Spring');
Map.addLayer(rgb_summer.mean(), visualization, 'RGB Summer',false);
var SensitivityBandVis = {
  min: 0,
  max: 50,
  palette: ['blue', 'yellow', 'green']
};
//Map.addLayer(rgb_spring.mean(), visualization, 'RGB Spring');
//Map.addLayer(sensitivity, SensitivityBandVis, 'Sensitivity');

var MvBandVis = {
  min: 0,
  max: 1,
  palette: ['blue', 'yellow', 'green']
};

var DeltaBandVis = {
  min: -1,
  max: 1,
  palette: ['blue', 'yellow', 'green']
};
Map.addLayer(Mv_primavera, MvBandVis, 'Mv Spring',false);
Map.addLayer(Mv_estate, MvBandVis, 'Mv Summer',false);
Map.addLayer(delta, DeltaBandVis, 'Mv Delta',true,0.5);

martedì 15 novembre 2022

Soil Moisture da Sentinel 1

Frugando nei custom scripts di Sentinel 1 ho trovato questo algoritmo che stima l'umidita' del suolo da dati radar in base all'intensita' di  backscattering. Si tratta di una analisi multitemporale prendendo il valore massimo e minimo di una serie di almeno 3 anni; questo approccio permette di ridurre di effetti di scabrosita' del terreno (si presuppone che nel tempo le condizioni nel pixel si mantengano stabili con la sola variazione dell'umidita', anche l'eventuale crescita di vegetazione influenza i valori di backscattering)

Il principio fisico su cui si basa l'algoritmo e' la modifica della costante dielettrica del terreno sulla base del contenuto di acqua (bibliografia)

In alcuni casi e' possibile identificare una legge di correlazione tra intensita' di backscattering ed umidita' tramite verit' a terra ma l'algoritmo del custom script definisce solo categorie di variazioni e non valori assoluti di umidita', Alcuni autori hanno ritrovato valori di correlazione fino all'80% ma in bibliografia vengono indicati casi di correlazione del tutto assente

Esempio di correlazione con verita' a terra

Per testare l'algoritmo su Sentinel Playground si puo' andare a questo link

Ho provato a metterlo alla prova nell'area del Mugello perche' nel settore Nord nei terrazzi alluvionali e' frequente vedere, anche grazie ai lavori agricoli che rimuovono la vegetazione, zone umide anche in piena estate in corrispondenza di passaggi di litologia (sono frequenti anche sorgenti di contatto lungo le incisioni dei torrenti)




Volevo vedere se i dati radar avevano una corrispondenza con la geologia utilizzando la cartografia CARG che al momento ha il miglior dettaglio nella zona di interesse



I risultati sono nelle due immagini rappresentativi del periodo secco (agosto) e quello piu' umido (novembre). E sono un po' perplesso. I dati non sono casuali perche' comunque c'e' coerenza tra i due periodi. Le aree blu (umide) si attestano non alla base degli impluvi ma sulle creste dei terrazzi. Si vede un modesto controllo litologico sulla distribuzione dell'umidita' ma quello piu' evidente dall'elaborazione e' quello di tipo topografico

Ho pensato anche ad un errore nella georefenziazione delle immagini di Sentinel Playground (che ho dovuto fare a mano perche' non mi ha mai funzionato il download automatico) ma dopo vari controlli lo escluderei

Il progetto QGIS e' scaricabile qui


[1] Wagner, W., Lemoine, G., Borgeaud, M. and Rott, H., 1999. A study of vegetation cover effects on ERS scatterometer data. IEEE Transactions on Geoscience and Remote Sensing, 37(2), pp.938-948.

[2] B. Bauer-Marschallinger et al., "Toward Global Soil Moisture Monitoring With Sentinel-1: Harnessing Assets and Overcoming Obstacles," in IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 1, pp. 520-539, Jan. 2019.

Tour de France Firenze Gran Depart

  Poagcar Wout Van Aert Vingegaard       Van Der Poel