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

Interferogramma Sentinel 1

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

https://asf.alaska.edu/how-to/data-recipes/create-an-interferogram-using-esas-sentinel-1-toolbox/

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 

 

Aggiornamento automatico temi Geoserver

Questo post serve a dare una procedura per effettuare un aggiornamento automatico di un tematismo raster con Geoserver

Prima di tutto si utilizza per semplicita' un container docker per lanciare il server

sudo docker run -it -p 8085:8080 --mount src="/home/luca/geoserver_data",target=/opt/geoserver_data/,type=bind docker.osgeo.org/geoserver:2.21.1

il server viene esposto su porta 8085 ed i dati sono da mettere nel folder /home/luca/geoserve

Per amministrare Geoserver si punta questo indirizzo

http://localhost:8085/geoserver/web

si aggiunge uno Store selezionando Geotiff, si seleziona lo WorkSpace (se non esiste va configurato prima), in Connection Parameters Url si clicca Browse e si seleziona il file raster di interesse

Dopo di cio si deve cliccare Enabled e Publish per rendere effettiva la visualizzazione

Per visualizzare il layer WMS Geoserver in Qgis si aggiunge nuovo layer WMS con indirizzo 

http://localhost:8085/geoserver/wms/

e poi si seleziona il tema da visualizzare

  

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>

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

Compressione immagine per Aruco Tags

Per il progetto Aruco ho provato ad acquisire un fotogramma da una webcam  e salvarlo in diversi formati con diverso grado di compressione. Di sotto si riportano i risultati (ingrandire le immagini per vedere gli artefatti di compressione)

import cv2
camera = cv2.VideoCapture(0)
camera.set(cv2.CAP_PROP_FRAME_WIDTH,1280)
camera.set(cv2.CAP_PROP_FRAME_HEIGHT,720)
compression_params = [cv2.IMWRITE_PNG_COMPRESSION, 0] 
while True:
    return_value,image = camera.read()
    gray = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)
    cv2.imshow('image',gray)
    if cv2.waitKey(1)& 0xFF == ord('s'):
        cv2.imwrite('test_0.png',image,[cv2.IMWRITE_PNG_COMPRESSION,0])
        cv2.imwrite('test_3.png',image,[cv2.IMWRITE_PNG_COMPRESSION,3])
        cv2.imwrite('test_9.png',image,[cv2.IMWRITE_PNG_COMPRESSION,9])
        cv2.imwrite('test_100.jpg',image,[cv2.IMWRITE_JPEG_QUALITY,100])
        cv2.imwrite('test_060.jpg',image,[cv2.IMWRITE_JPEG_QUALITY,60])
        cv2.imwrite('test_060.ppm',image,[cv2.IMWRITE_PXM_BINARY,1])
        break
camera.release()
cv2.destroyAllWindows()

PNG compressione 0

PNG compressione 3

PNG compressione 9

JPG quality 60

JPG quality 100

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

Geotiff da assets in EarthEngine

In Earth Engine e' possibile visualizzare dei Geotiff caricandoli in assets 

I geotiff devono essere puri e non tiff con allegato il world file in formato tfw.

Per la conversione da tfw a geotiff si puo' usare Gdal Translate includendo anche il sistema di riferimento

(quello che segue e' un esempio di conversione di un foglio CARG Carta geologica regionale Toscana)

gdal_translate -a_srs EPSG:3003 -of GTiff 264020.tif 264020_.tif

Le immagini che sono in assets non sono immediatamente leggibili dalla versione app del progetto.

Devono essere esplicitamente impostati i permessi di lettura nelle proprieta' dell'asset

Per visualizzare file prendendolo dagli assets si usa la seguente sintassi

var geologia = ee.Image('projects/ee-lucainnoc/assets/264020');

Map.addLayer(geologia, {},'Geologia',false,0.5);

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

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

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.

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

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

Truecolor Sentinel2 con RasterIO

Truecolor da dati Sentinel 2 usando RasterIO 

L'immagine e' stata corretta nel contrasto per renderla piu' leggibile

import rasterio

from rasterio import plot

imagePath = './immagini/'

band2 = rasterio.open(imagePath+'B2.jp2', driver='JP2OpenJPEG') #blue

band3 = rasterio.open(imagePath+'B3.jp2', driver='JP2OpenJPEG') #green

band4 = rasterio.open(imagePath+'B4.jp2', driver='JP2OpenJPEG') #red

band2s = band2.read(1)

band3s = band3.read(1)

band4s = band4.read(1)

bandg2=2.5*(band2s.astype(float))

bandg3=2.5*(band3s.astype(float))

bandg4=2.5*(band4s.astype(float))

trueColor = rasterio.open('./immagini/SentinelTrueColor2.tiff','w',driver='Gtiff',

width=band4.width, height=band4.height,

count=3,

crs=band4.crs,

transform=band4.transform,

dtype=band4.dtypes[0]

)

trueColor.write(bandg2,3) #blue

trueColor.write(bandg3,2) #green

trueColor.write(bandg4,1) #red

trueColor.close()

Sentinel 2 da SciHub

Nel precedente post avevo provato a scaricare i dati Sentinel da un bucket Google Cloud

In realta' ho trovato una procedura piu' lineare usando OpenHub di https://scihub.copernicus.eu/ 

E' necessario un account ma non e' legato a nessuna carta di credito ed il download e' libero

Per effettuare la selezione ed il download si usa lo script dhusget (Altre istruzioni a questo link )

una stringa di richiesta e' del tipo

./dhusget_0.3.8.sh -u username -p password -m Sentinel-2 -T S2MSI2A -t96 -o product -c 10.6,43.5:11.3,44 -n 1 -w 5 -O ./immagini

in cui viene effettuata la richeista delle immagini che contengano Firenze (vedi coordinate) acquisite nelle ultime 96 ore della missione Sentinel 2

Si possono specificare anche intervalli temporali tramite intervalli di date (in questo caso tra il 10/8/2022 ed il 12/8/2022

./dhusget_0.3.8.sh -u c1p81 -p xxxxxx@xxxxx -m Sentinel-2 -T S2MSI2A -S 2022-08-10T00:00:00.000Z -E 2022-08-12T00:00:00.000Z -n 1 -w 5 -c 10.6,43.5:11.3,44 -o product  -O ./immagini

Si possono scaricare sia il livello L1C che L2A .. di seguito si riporta il contenuto del file zip

I file zip vanno da circa 800 Mb ad oltre 1.1 Gb

Per selezionare solo un prodotto si puo' usare lo switch 

-T S2MSI1C per L21C

-T S2MSI2A per L2A

Ho notato che il sistema spesso tronca la connessione e non termina in modo corretto il download...se si abilita la funzione di checksum al termine non funziona mai..disabilitandolo invece si ottengono dei download corretti

LEVEL 1C

tutte le bande ognuna alla sua risoluzione nativa + truecolor a 10 m

LEVEL 2A 
risoluzione 10 m
AOT : Aerosol Optical Thickness
WVP : Water vapor 

TCI : Truecolor

Bande B2,B3,B4,B8

risoluzione 20 m
AOT : Aerosol Optical Thickness
WVP : Water vapor 

TCI : Truecolor

SCL : Scene classification map

Bande B1,B2,B3,B4,B5,B6,B7,B11,B12,B8A

risoluzione 60 m
AOT : Aerosol Optical Thickness
WVP : Water vapor 

TCI : Truecolor

SCL : Scene classification map

Bande B1,B2,B3,B4,B5,B6,B7,B9,B11,B12,B8A

con il medesimo sistema si possono ottenere i dati di Sentinel 1 SLC, OCN, GRDH e RAW

Sentinel 2 in Google Cloud

 Le immagini Sentinel2 sono rilasciate in modo gratuito ma i repository delle immagini le forniscono con una interfaccia utente interattiva mentre e' a pagamento l'accesso automatico via API

In modo abbastanza casuale mi sono imbattuto in questo bucket su Google Cloud

https://console.cloud.google.com/storage/browser/gcp-public-data-sentinel-2/tiles/33/U/UP;tab=objects?prefix=&forceOnObjectsSortingFiltering=false

che puo' essere utilizzato mediante questa libreria Python (per la selezione e download)

https://github.com/QuantuMobileSoftware/sentinel2tools

per l'installazione si usa

pip3 install 
git+https://github.com/QuantuMobileSoftware/sentinel2tools.git@ca232cb66106db6cac729defdab91aad9aecb15b.

Lo scripthe puo' essere utilizzato mediante questa libreria Python (per la selezione e download)

Il criterio di selezione geografica avviene tramite un geojson che puo' essere creato mediante il servizio on line geojson.io

from sentinel2download.downloader import Sentinel2Downloader

from sentinel2download.overlap import Sentinel2Overlap

if __name__ == '__main__':

verbose = True

aoi_path = "./firenze_geojson_utm.json"

overlap = Sentinel2Overlap(aoi_path, verbose=verbose)

tiles = overlap.overlap()

print(f"Overlapped tiles: {tiles}")

api_key = f"./api_key.json"

loader = Sentinel2Downloader(api_key, verbose=verbose)

product_type = 'L2A' # or L1C

start_date = "2022-09-01"

end_date = "2022-09-30"

output_dir = './utm'

cores = 3

BANDS = {'TCI', 'B04', }

CONSTRAINTS = {'NODATA_PIXEL_PERCENTAGE': 15.0, 'CLOUDY_PIXEL_PERCENTAGE': 10.0, }

loaded = loader.download(product_type,

tiles,

start_date=start_date,

end_date=end_date,

output_dir=output_dir,

cores=cores,

bands=BANDS,

constraints=CONSTRAINTS)

print(f"Load information")

for item in loaded:

print(item)

print("Execution ended")

Visto che i dati sono in un bucket di Gooogle Cloud e' richiesta anche 'autenticazione a GCloud. 

Per avere il file json di autenticazione si deve aprire la consolle di GCloud, creare un progetto e seguire le istruzioni al link seguente

https://cloud.google.com/api-keys/docs/get-started-api-keys

Le immagini vengono scaricate in formato JP2, una per banda, a 16 bit e georiferite in UTM (attenzione a questo dato). Oltre alle 12 bande sono presenti un truecolor (TCI) 

Viene scaricata tutta l'immagine che comprende la selezione geografica. Per effettuare il taglio si possono usare le GDAL. Le coordinate devono essere UTM

from osgeo import gdal

upper_left_x = 674882

upper_left_y = 4854678

lower_right_x = 688882

lower_right_y = 4844623

window = (upper_left_x,upper_left_y,lower_right_x,lower_right_y)

gdal.Translate('output_tci.jp2', 'tci.jp2', projWin = window)