GRIB Arpege MeteoFrance

Un'altra fonte di download di previsioni meteo specifico per l'area Europa e' il servizio Arpege con Grib files che si possono scarica da 

https://donneespubliques.meteofrance.fr/?fond=produit&id_produit=130&id_rubrique=51

Le previsioni sono orarie fino a 114 ore con una risoluzione di 0.1 gradi

dal sito non e' possibile selezionare la finestra geografica dei dati. Per effettuare una selezione sul file GRIB si puo' usare il parametro -lon ...il problema e' che questo switch non e' compatibile con lo switch -cvs per cui per fare l'estrazione si deve usare il piping

di seguito uno script che scarica in automatico i file disponibili, estra il valore temperatura per la coordinata selezionata e crea un file csv 

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

#!/bin/bash
rm *.grib > /dev/null 2>&1
rm *.csv > /dev/null 2>&1

today=`date +%Y-%m-%d`
today2=`date +%Y%m%d`

lon="11.8"
lat="46.4"
 
wget -O 00_12.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=00H12H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 13_24.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=13H24H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 25_36.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=25H36H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 37_48.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=37H48H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 49_60.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=49H60H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 61_72.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=61H72H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 73_84.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=73H84H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 85_96.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=85H96H&referencetime=${today}T00:00:00Z&format=grib2" 
wget -O 97_102.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=97H102H&referencetime=${today}T00:00:00Z&format=grib2" 
#wget -O 103_114.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZp6y0leVHU__&model=ARPEGE&grid=0.1&package=SP1&time=103H114H&referencetime=${today}T00:00:00Z&format=grib2" 
.

/wgrib2 -for 102:114 -s -lon ${lon} ${lat} 00_12.grib2 > dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 13_24.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 25_36.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 37_48.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 49_60.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 61_72.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 73_84.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 85_96.grib2 >> dati.csv
./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 97_102.grib2 >> dati.csv
#./wgrib2 -for 102:114 -s -lon ${lon} ${lat} 103_114.grib2 >> dati.csv

cut -d, -f1 --complement dati.csv > dati2.csv

cut -d, -f1 --complement dati2.csv > temp.csv

sed -i -e 's/val=//g' temp.csv

GFS via NOMAD

 Se non si ha necessita' di scaricare tutto il dataset di GFS si puo' effettuare una query di selezione utilizzando la URL sottostante

https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?dir=%2Fgfs.20220724%2F00%2Fatmos

dopo si puo' scegliere il campo desiderato. Verra' scaricato un file GRIB 

il tracciato record completo e' https://www.nco.ncep.noaa.gov/pmb/products/gfs/gfs.t00z.pgrb2.0p25.f000.shtml

Il download puo' essere automatizzato fornendo i parametri sulla URL (in questo caso si richiede il parametro TMP 2 m above ground nella finestra latitudine 41-46, longitudine 10-12 per la previsione F000 del giorno 24/07/2022

https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f000&lev_2_m_above_ground=on&var_TMP=on&subregion=&leftlon=10&rightlon=12&toplat=46&bottomlat=41&dir=%2Fgfs.20220724%2F00%2Fatmos

Attenzione che EarthEngine di Google esprime i valori in C mentre da questa fonte i dati sono in K

se interessa la precipitazione totale si puo' usare il parametro PRATE (Precipitation is the total liquid content of rain, snow, sleet (ice pellets), and freezing rain which falls out of the atmosphere in a given time period. Precipitation rate is the rate at which precipitation is falling at a given momen) espresso in kg/m2/sec (1Kg/m2 corrisponde ad 1l/m2 quindi 1mm/m2, per avere i dati orari in mm/ora basta moltiplicare per 3600). Anche in questo caso vi sono differenze tra come espressi i dati (Kg/m2/s contro Kg/m2)

https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.anl&all_lev=on&var_PRATE=on&leftlon=7.75&rightlon=8&toplat=45&bottomlat=44.5&dir=%2Fgfs.20220724%2F00%2Fatmos

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

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