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

mercoledì 31 maggio 2023

Batimetria Venezia da Landsat con Decision Tree

Visto lo scarso successo del metodo analitico visto qui ho provato un metodo alternativo citato spesso il letteratura che prevede la determinazione della batimetria da dati ottici usando la rete neurale decision tree

Per cambiare ho preso i dati Landsat sulla Laguna di Venezia (la batimetria della Laguna di Venezia e' molto piu' variabile rispetto a quella della lagunan di Orbetello) usando come dato di verita' a terra questo geotiff http://cigno.ve.ismar.cnr.it/layers/geonode%3Alag02_okart_export

Su Qgis sono stati importate le bande RGB,NIR di Landsat8  e dopo avere creato una griglia regolare di punti spaziati di 30 m e' stata estratta una tabella con i valori di riflettanza per ogni banda e la profondita'


 Una volta ottenuto il file CSV i dati di batimetria sono stati divisi in classi di 0.5 m per rendere possibile la successiva elaborazione tramite un semplice script in GO 


package main

import (
    "encoding/csv"
    "fmt"
    "io"
    "log"
    "math"
    "os"
    "strconv"
)

func main() {
    f, err := os.Open("tree.csv")
    if err != nil {
        log.Fatal(err)
    }
    defer f.Close()
    csvReader := csv.NewReader(f)
    for {
        rec, err := csvReader.Read()
        if err == io.EOF {
            break
        }
        if err != nil {
            log.Fatal(err)
        }
        tt, _ := strconv.ParseFloat(rec[0], 'f')
        t := math.Abs(tt)
        fmt.Print(math.Floor(t / 0.5))
        fmt.Println("," + rec[1] + "," + rec[2] + "," + rec[3] + "," + rec[4])

    }
}


il file finale ha un formato del tipo

 classe,red,nir,gren,blue
3,6707,5566,8241,9397
3,6714,5575,8221,9375
3,6696,5573,8184,9369
3,6665,5577,8144,9331
3,6638,5584,8089,9287
3,6636,5568,8080,9281

In totale sono emerse 23 classi per un totale di 16000 punti


Come si vede la distribuzione e' fondamentalmentre asimmetrica con le prima 4 classi che rappresentano la gran parte dei dati. E' stato effettuato un taglio dei dati alle prime 6 classi

 Il file CSV e' stato utilizzato  all'interno dello script seguente in Colab (si tratta di un semplice adattamento del programma di esempio di Tensorflow per la rete  decision tree)


# -*- coding: utf-8 -*-
"""decisiontree_venezia.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1f8HLdjHiDNvpplKzNOVrXuIXrKTay0TS
"""

!pip install tensorflow_decision_forests

import numpy as np
import pandas as pd
import tensorflow_decision_forests as tfdf

path = "/content/classi_nobat_ridotto.csv"
pandas_dataset = pd.read_csv(path)

# Display the first 3 examples.
pandas_dataset.head(3)

np.random.seed(1)
# Use the ~10% of the examples as the testing set
# and the remaining ~90% of the examples as the training set.
test_indices = np.random.rand(len(pandas_dataset)) < 0.1
pandas_train_dataset = pandas_dataset[~test_indices]
pandas_test_dataset = pandas_dataset[test_indices]

print("Training examples: ", len(pandas_train_dataset))

print("Testing examples: ", len(pandas_test_dataset))

tf_train_dataset = tfdf.keras.pd_dataframe_to_tf_dataset(pandas_train_dataset, label='classe')
model = tfdf.keras.CartModel()
model.fit(tf_train_dataset)

tfdf.model_plotter.plot_model_in_colab(model, max_depth=10)

model.compile("accuracy")
print("Train evaluation: ", model.evaluate(tf_train_dataset, return_dict=True))

tf_test_dataset = tfdf.keras.pd_dataframe_to_tf_dataset(pandas_test_dataset, label='classe')
print("Test evaluation: ", model.evaluate(tf_test_dataset, return_dict=True))

from keras.utils import plot_model

plot_model(
    model,
    to_file='model.png',
    show_shapes=False,
    show_dtype=False,
    show_layer_names=True,
    rankdir='TB',
    expand_nested=False,
    dpi=96,
    layer_range=None,
    show_layer_activations=False,
    show_trainable=False
)

!pip install keras-tuner

import keras_tuner as kt

def build_model(hp):
  model = tfdf.keras.CartModel(
      min_examples=hp.Choice("min_examples",
          # Try four possible values for "min_examples" hyperparameter.
          # min_examples=10 would limit the growth of the decision tree,
          # while min_examples=1 would lead to deeper decision trees.
         [1, 2, 5, 10]),
      validation_ratio=hp.Choice("validation_ratio",
         # Three possible values for the "validation_ratio" hyperparameter.
         [0.0, 0.05, 0.10]),
      )
  model.compile("accuracy")
  return model

tuner = kt.RandomSearch(
    build_model,
    objective="val_accuracy",
    max_trials=10,
    directory="/tmp/tuner",
    project_name="tune_cart")

tuner.search(x=tf_train_dataset, validation_data=tf_test_dataset)
best_model = tuner.get_best_models()[0]

print("Best hyperparameters: ", tuner.get_best_hyperparameters()[0].values)
# >> Best hyperparameters:  {'min_examples': 2, 'validation_ratio': 0.0}

model = tfdf.keras.CartModel(min_examples=2, validation_ratio=0.0)
model.fit(tf_train_dataset)

model.compile("accuracy")
print("Test evaluation: ", model.evaluate(tf_test_dataset, return_dict=True))

tfdf.model_plotter.plot_model_in_colab(model, max_depth=10)

Training examples: 13463 Testing examples: 1520


 



al termine si ha una accuracy della porzione di dati di test superiore al 70%

14/14 [==============================] - 1s 8ms/step - loss: 0.0000e+00 - accuracy: 0.7761 Train evaluation: {'loss': 0.0, 'accuracy': 0.7761271595954895} 2/2 [==============================] - 0s 19ms/step - loss: 0.0000e+00 - accuracy: 0.7388 Test evaluation: {'loss': 0.0, 'accuracy': 0.7388157844543457} 



In generale la rete nel caso di incertezza di una classe seleziona sempre classi contingue

 

mercoledì 17 maggio 2023

Pansharpening Landsat con Qgis

 

Si caricano i raster delle bande 2,3,4 (BGR) e la banda 8 pancromatica

Si apre il menu Raster/Miscellanea/Crea Raster Virtuale selezionando come input le bande ottiche e spuntando il flag Place each input file in a separate band

A questo punto dal toolbox si seleziona pansharpening indicando come Dataset spettrale il raster virtuale creato in precedenza e come dataset pancromatico la banda 8




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


martedì 23 agosto 2022

Landsat 8 truecolor gif con Google Earth Engine

 Una prova di LandSat 8 truecolor con Google Earth Engine...stavo cercando i periodi con neve su questo settore dell'Appennino ma non e' mi e' stato molto di aiuto (il periodo va da agosto 2013 a agosto 2022....28 immagini con lista completa al termine del post)


var geometry = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.11590875374765, 44.54755512745409],
          [10.11590875374765, 44.39388022912829],
          [10.329541337365814, 44.39388022912829],
          [10.329541337365814, 44.54755512745409]]], null, false);
var visualization = {"bands":["B4","B3","B2"],"min":0,"max":0.3};
// Applies scaling factors.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}
// Landat 8 surface reflection data
var L8Coll = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(geometry).sort('DATE_ACQUIRED')
    .map(applyScaleFactors)
    .filterMetadata('CLOUD_COVER', 'less_than', 1)
    .map(function(image){return image.clip(geometry)});
    
print(L8Coll)
var L8Coll=L8Coll.select(['SR_B2', 'SR_B3', 'SR_B4','SR_B5', 'SR_B6', 'SR_B7','ST_B10'] , ['B2', 'B3', 'B4','B5', 'B6', 'B7','B10'])
print(L8Coll)
Map.addLayer(L8Coll.first(), visualization, 'L8 2013 ');   
Map.addLayer(L8Coll.sort('CLOUD_COVER').first(), visualization, 'L8 Best ');   
Map.centerObject(L8Coll)
print('L8Coll Least Cloud Cover',L8Coll.sort('CLOUD_COVER').first())
print('Landsat Collection Date Acquired',L8Coll.aggregate_array('DATE_ACQUIRED'))
print('Landsat Collection Cloud Cover',L8Coll.aggregate_array('CLOUD_COVER'))
print('Landsat Collection Path/Row',L8Coll.aggregate_array('WRS_PATH'),
      L8Coll.aggregate_array('WRS_ROW'))
var L8LeastCC=L8Coll.sort('CLOUD_COVER').first()
print('L8 Least Cloud Cover',L8LeastCC)
print(ee.Date(L8Coll.first().get('system:time_start')).format("yyyy-MM-dd"))
//+++++++++++++++++++++
/*
// Load a Landsat 8 collection for a single path-row.
var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filter(ee.Filter.eq('WRS_PATH', 179))
  .filter(ee.Filter.eq('WRS_ROW', 34));
*/
// This function adds a band representing the image timestamp.
var addTime = function(image) {
  return image.addBands(image.metadata('system:time_start'));
};
// Map the function over the collection and display the result.
var L8CollT=L8Coll.map(addTime)
print('L8coll system Time',L8CollT);

print(ee.Date(L8Coll.first().get('system:time_start')).format("yyyy-MM-dd"))
var doy = ee.Date(L8Coll.first().get('system:time_start')).getRelative('day', 'year');
print('doy',doy)
//++++++++++++++++++++++

// Create RGB visualization images for use as animation frames.
var rgbVis = L8Coll.map(function(img) {
  return img.visualize(visualization).clip(geometry)
});
print('rgbVis',rgbVis)
// Define GIF visualization parameters.
var gifParams = {
  'region': geometry,
  'dimensions': 800,
  'crs': 'EPSG:3857',
  'framesPerSecond': 1//,  'format': 'gif'
};
// Print the GIF URL to the console.
print(rgbVis.getVideoThumbURL(gifParams));
// Render the GIF animation in the console.
print(ui.Thumbnail(rgbVis, gifParams));
////+++++++++++++++

0:
2013-08-03
1:
2013-08-19
2:
2013-09-04
3:
2014-06-10
4:
2016-07-17
5:
2016-08-27
6:
2017-03-14
7:
2017-04-08
8:
2017-06-11
9:
2017-08-05
10:
2017-08-30
11:
2018-06-30
12:
2019-01-15
13:
2019-06-24
14:
2019-07-26
15:
2019-09-12
16:
2019-09-21
17:
2020-04-07
18:
2020-04-23
19:
2020-07-21
20:
2020-08-22
21:
2021-03-02
22:
2021-10-28
23:
2022-02-08
24:
2022-03-28
25:
2022-04-29
26:
2022-07-02
27:
2022-08-03

Geologi

  E so anche espatriare senza praticamente toccare strada asfaltata