DHOB (IU5SGN)

venerdì 22 luglio 2022

EarthEngine con Python

Una delle limitazione di Earth Engine Editor e' quello di salvare i dati solo su cloud. Tale limitazione puo' essere aggirata mediante Geemap che permette di scarica i dati come array numpy o come geotiff

Prima di poter procedere con Geemap si deve configurare l'autenticazione con EarthEngine mediante

pip install earthengine-api --upgrade

=====================================================

import ee
import geemap
import numpy as np
import os

ee.Authenticate()
ee.Initialize()
img = ee.Image('NOAA/GFS0P25/2022072100F002').select(['temperature_2m_above_ground'])
aoi = ee.Geometry.Polygon([[[11.0, 43.0], [11.0, 43.25], [11.25, 43.25], [11.25, 43.0]]], None, False)

# converte una immagine in un array numpy
rgb_img = geemap.ee_to_numpy(img, region=aoi)
print(rgb_img.shape)
print(rgb_img)

# scarica una image collection in una serie di Geotiff
aoi = ee.Geometry.Point(11.0, 43.0)
region = aoi.buffer(1000)
collection = (
ee.ImageCollection('NOAA/GFS0P25')
.filterBounds(aoi)
.filterDate("2022-07-21")
.select('temperature_2m_above_ground')
)
out_dir='/home/luca/'
geemap.ee_export_image_collection(collection, out_dir=out_dir, region=region)

=====================================================

per scaricare una intera immagine si puo' utilizzare la libreria Geedim (attenzione le versioni precedenti alla 1.3 hanno problemi per il download di grandi dimensioni)

un esempio puo' essere il seguente per ricercare le immagini disponibili

geedim search -c NOAA/GFS0P25 -s 2022-07-21 -e 2022-07-22 --bbox 43.0 11.0 44.0 12.0

il bbox indica l'angolo in basso/sinistra alto/destra WGS84

per scaricare l'immagine si puo' usare la sintassi sottostante

geedim download --id NOAA/GFS0P25/2022072118F000 --crs EPSG:4326 --bbox 43.0 11.0 44.0 12.0 --scale 1000 -o

-o : overwrite il file se esiste

il geotiffi contiene 7 bande (1: temperature 2 m above ground, 2: umidita' specifica, 3: umidita' relativa, 4,5 : vento, 6: acqua precipitabile in atmosfera, 7: maschera)

NOAA GFS correlazione temperatura con dati a terra

Una prova speditiva di correlazione tra dati a terra e dati di previsione dal servizio NOAA GFS

E' stata selezionata (per motivi non casuali) la stazione meteo di Punta Rocca con tutti i problemi di correlare un dato puntuale della stazione meteo con il dato esteso del pixel del dato da satellite e l'orografia dell'area

Prendendo i dati in modo acritico si vede che in valore assoluto non c'e' corrispondenza tra satellite e terra ma effettuando uno scatterplot si vede una correlazione non trascurabile (considerando anche l'incertezza di un dato di previsione con una misura reale)

I dati sono dati orari (e' stato corretto l'orario UTC del satellite rispetto all'orario solare italiano della stazione meteo) e sono scalati dal 15 al 22 luglio 2022

giovedì 21 luglio 2022

Google Earth Cloud Engine NOAA GFS

Invece di procedere in modo tradizionale come indicato nel post precedente si possono interrogare i dati GFS NOAA anche da Google Earth Engine con il vantaggi che i dati sono gia' scaricati e pronti all'elaborazione (con lo svantaggio che i risultati possono essere esportati in automatico solo su cloud e non sulla macchina locale)

=================================================

var palettes = require('users/gena/packages:palettes');
var palette = palettes.colorbrewer.RdYlGn[9];

var region = ee.Geometry(Map.getBounds(true))
// Strip time off current date/time
var today = ee.Date(new Date().toISOString().split('T')[0])
var image;
var collection = ee.ImageCollection('NOAA/GFS0P25')
.select('temperature_2m_above_ground')
.filterDate(today.advance(-1, 'day'), today.advance(1, 'day'))
.sort('system:time_start', false)
console.log(collection);
// per i successivi 3 giorni le previsioni sono orarie
// le previsioni vengono ricalcolate ogni 6 ore
var forecasts = ee.ImageCollection(
[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23, 24,
25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,
49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,63,64,65,66,67,68,69,70,71].map(forecast)
//[12, 24, 36,48,60, 72].map(forecast)
)
Map.setCenter(11.8537, 46.4297, 12);
Map.setOptions('SATELLITE');
Map.setControlVisibility(false, false, true, true, true, false, false)
Map.addLayer(forecasts, {min: 0, max: 30,palette: palette,opacity: 0.5}, 'Temperature forecast', true)
var textStyle = {
color: 'grey',
fontName: 'arial',
fontSize: 10,
bold: true,
italic: false
}
var chart = ui.Chart.image.series({
imageCollection: forecasts,
region: geometry,
reducer: ee.Reducer.mean(),
scale: forecasts.first().projection().nominalScale(),
xProperty: 'date',
}).setChartType('LineChart').setOptions({title: 'Forecast 72h Temperature 2m above ground (C)',hAxis: {title: 'Date',textStyle: textStyle},vAxis: {title: 'Temperature'}});

//print(chart)
var panel_chart = ui.Panel({style: {width: '400px'}})
.add(ui.Label('NOAA/GFS0P25'))
.add(chart);
ui.root.add(panel_chart);

function forecast(hours) {
image = collection
.filterMetadata('forecast_hours', 'equals', hours)
// Since colleciton is sorted descending, if there are more than
// one forecast, first() will give the latest.
.first()
var date = image.date().advance(hours, 'hours')
return image
.set('date', date.format())
}
var panel = ui.Panel();
panel.style().set({
width: '400px',
position: 'bottom-right'
});
//Map.add(panel);
Map.onClick(function(coords) {
//panel_chart.clear();
console.log(collection);
//Map.add(panel);
var location = 'lon: ' + coords.lon.toFixed(4) + ' ' +
'lat: ' + coords.lat.toFixed(4);
var click_point = ee.Geometry.Point(coords.lon, coords.lat);
var valore = image.reduceRegion(ee.Reducer.first(), click_point, 90).evaluate(function(val){
//console.log(valore);
//var demText = 'Habitat suitability: ' + val.b1;
//toolPanel.widgets().set(2, ui.Label(demText));
});
panel.widgets().set(1, ui.Label(location));
// Edit: To be temporary, the "loading..." panel number has to be the same as the demText panel number (changed from 1 to 2).
//panel.widgets().set(2, ui.Label("loading..."));
Map.layers().set(1, ui.Map.Layer(click_point, {color: 'FF0000'}));
});

Utilizzo del servizio NOAA GFS

I files vengono distribuiti in originale dal NOAA via FTP in formato Grib2 dall'indirizzo

ftp.ncep.noaa.gov

nel folder

/pub/data/nccf/com/gfs/v16.2/gfs.20220721/00/atmos

dove 20020721 e' la data che deve essere modificata

ed il nome del file e' nel formato

gfs.t00z.pgrb2.0p25.f000

la parte finale f000 varia da f000 a f384 ed indica l'ora di previsione

(ogni giorni vengono effettuati 4 run di calcolo alle 00,06,012,18 e vengono calcolate le previsioni per le successive 384 ore). Gli orari sono in UTC

si puo' eseguire in download anche via curl

data=$(date +%Y%m%d)
curl -Os ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/v16.2/gfs.$data/00/atmos/gfs.t00z.pgrb2.0p25.f000

c'e' da fare attenzione che ogni file e' di dimensioni oltre i 450 Mb

i campi dati sono descritti qui e sono scalati su una griglia di 0.25 gradi

i dati sono in formato GRIB2 e possono essere decodificati usando il programma WGRIB2 scaricabile da questo indirizzo https://ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/wgrib2.tgz

Per compilare su Linux si spacchetta il file tar e poi da linea di comando si digita

export CC=gcc
export FC=gfortran
export COMP_SYS=gnu_linux
make

per estrarre i dati dal file (per esempio il campo Temperature 2 m Above Ground e' il 580) nel formato csv

./wgrib2 ../../gfs.t00z.pgrb2.0p25.f000 -s -d 580 -csv dati.csv

per estrarre i dati ad una determinata posizione geografica

/wgrib2 ../../gfs.t00z.pgrb2.0p25.f000 -s -d 580 -lon 46.5 11.75

che riporta come risultato

580:418584888:d=2022072100:TMP:2 m above ground:anl::lon=46.500000,lat=11.750000,val=301.51

Un ultimo comodo comando puo' essere quello di trasformare il GRIB in un Geotiff mediante

gdal_translate -b 580 -a_srs EPSG:4326 gfs.t00z.pgrb2.0p25.f000 -of Gtiff test.geotiff

il parametro -b indica la banda che del file grib che si vuole esportare