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.

lunedì 14 novembre 2022

Mascherare le nuvole in Sentinel 2

Nei dati Sentinel 2 e' contenuta anche in QI_DATA/MSK_CLDPRB_20m.jp2 una maschera di probabilita' che il pixel sia relativo ad una nuova (con valore di probabilita' crescente)




Per fare una maschera ho provato a settare i valori maggiori di 50 a zero ed i valori inferiori a 1 per poi moltiplicare i valori della banda



import numpy as np
import rasterio
from rasterio import plot

band4 = rasterio.open('b4.jp2', driver='JP2OpenJPEG')
mask = rasterio.open('mask.jp2', driver='JP2OpenJPEG')

band4_ = band4.read(1)
mask_ = mask.read(1)

mask_[mask_ <= 50] = 1
mask_[mask_ > 50] = 0


ris = np.zeros(band4.shape, dtype=rasterio.uint32)
ris =band4_*mask_

kwargs = band4.meta
kwargs.update(driver='GTiff',dtype=rasterio.uint16,count=1)

maskimage = rasterio.open('maschera0_1.tiff','w',**kwargs )
maskimage.write(mask_,1)
maskimage.close()

b4mask = rasterio.open('b4mask.tiff','w', **kwargs)
b4mask.write(ris,1)
b4mask.close()

Tagliare tutte le bande Sentinel in un colpo solo

 Puo' essere utile tagliare un granule di Sentinel sulla base di uno shape file


Per fare cio' si prepara lo shape (io ho usato lo stesso sistema di riferimento dell'immagine Sentinel... non so se rasterio funziona anche con sistema di riferimento misti) 

Tutte le bande sono nel folder immagini mentre il risultato sara' salvato nel folder rst

Nonostante sia Python e' uno script molto veloce


import os
import fiona
import rasterio
from rasterio.mask import mask

bandPath = './immagini/'
bandNames = os.listdir(bandPath)

aoiFile = fiona.open('./taglio.shp')
aoiGeom = [aoiFile[0]['geometry']]

for band in bandNames:
rasterPath = os.path.join(bandPath,band)
rasterBand = rasterio.open(rasterPath)
outImage, outTransform = mask(rasterBand, aoiGeom, crop=True)
outMeta = rasterBand.meta
outMeta.update({"driver": 'JP2OpenJPEG',
"height": outImage.shape[1],
"width": outImage.shape[2],
"transform": outTransform})
outPath = os.path.join('./rst',band)
outRaster = rasterio.open(outPath, "w", **outMeta)
outRaster.write(outImage)
outRaster.close()

sabato 12 novembre 2022

NDVI Sentinel 2 con RasterIO

il formato di output e' stato modificato in Geotiff perche' JP2 non permette di salvare in float32 


import numpy as np
import rasterio
from rasterio import plot

imagePath = './immagini/'
band4 = rasterio.open(imagePath+'B4.jp2', driver='JP2OpenJPEG')
band8 = rasterio.open(imagePath+'B8.jp2', driver='JP2OpenJPEG')

band4s = band4.read(1)
band8s = band8.read(1)

RED=band4s.astype(float)
NIR=band8s.astype(float)
np.seterr(divide = "ignore", invalid = "ignore")

ndvi = np.zeros(RED.shape, dtype=rasterio.float32)
ndvi = (NIR-RED) / (NIR+RED)


ndvi[ndvi > 1] = np.nan
ndvi[ndvi < -1] = np.nan

kwargs = band4.meta
kwargs.update(driver='GTiff',dtype=rasterio.float32,count=1)

ndviimage = rasterio.open('./immagini/SentinelNDVI2.tiff','w',**kwargs )
ndviimage.write(ndvi,1)
ndviimage.close()


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