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