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