Previsioni temperatura con GFS NOAA ed Arpege

Ho provato a mettere a confronto le previsioni meteo dei servizi NOAA GFS ed Arpege visti nei precedenti post per una localita' italiana confrontando i dati con dati di verita' a terra dati da stazione meteo. Per rendere piu' agevole l'analisi il download ed il pretrattamento dei dati in CSV e' stato automatizzato  

Script per NOAA GFS

-----------------------------------------------------------------------------------

dove="/home/luca/gfs/"
rm $dove*.grib > /dev/null 2>&1
rm $dove*.csv > /dev/null 2>&1
rm /var/www/html/luca/gfs/completo.csv > /dev/null 2>&1
mv $dove*.txt ${dove}backup/ > /dev/null 2>&1

today=`date +%Y%m%d`
echo $today
# 11.8537, 46.4297
llon="11.75"
llat="46.25"
ulon="11.99"
ulat="46.49"
for i in $(seq -f "%03g" 0 120)
do
  echo $i
  wget -q -O ${dove}f${i}.grib "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f$i&lev_2_m_above_ground=on&var_TMP=on&subregion=&leftlon=${llon}&rightlon=${ulon}&toplat=${ulat}&bottomlat=${llat}&dir=%2Fgfs.$today%2F00%2Fatmos"
done
for i in $(seq -f "%03g" 123 3 384) #123 126 129
do
  echo $i
  wget -q -O ${dove}f${i}.grib "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f$i&lev_2_m_above_ground=on&var_TMP=on&subregion=&leftlon=${llon}&rightlon=${ulon}&toplat=${ulat}&bottomlat=${llat}&dir=%2Fgfs.$today%2F00%2Fatmos"
done

for i in $(seq -f "%03g" 0 120)
do
    echo "Conversione grib${i}"
  ${dove}wgrib2 ${dove}f${i}.grib -csv ${dove}f${i}.csv
done
for i in $(seq -f "%03g" 123 3 384)
do
  echo "Conversione grib${i}"
  ${dove}wgrib2 ${dove}f${i}.grib -csv ${dove}f${i}.csv
done
for i in $(seq -f "%03g" 0 120)
do
  cat ${dove}f${i}.csv >> ${dove}${today}_completo.txt
done
for i in $(seq -f "%03g" 123 3 384)
do
  cat ${dove}f${i}.csv >> ${dove}${today}_completo.txt
done
#aggiunge l'intestazione al csv
sed '1 s/^/data1,data2,variabile,unita,lon,lat,valore\n/' ${dove}${today}_completo.txt > ${dove}${today}_completo_1.txt
cp ${dove}${today}_completo_1.txt /var/www/html/luca/gfs/completo.csv

---------------------------------------------------------------------------------

Script per Arpege (il servizio necessita di una API Key nel campo Token della URL ...ho modificato per non renderla utilizzabile la mia)

---------------------------------------------------------------------------------

lon="11.8"

lat="46.4"

 

wget -O ${dove}00_12.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=00H12H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}13_24.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=13H24H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}25_36.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=25H36H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}37_48.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=37H48H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}49_60.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=49H60H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}61_72.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=61H72H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}73_84.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBxxxxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=73H84H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}85_96.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxxxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=85H96H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}97_102.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxxxxxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=97H102H&referencetime=${today}T00:00:00Z&format=grib2" 

${dove}wgrib2 -for 102:114 -s -lon ${lon} ${lat} ${dove}00_12.grib2 > ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}13_24.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}25_36.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}37_48.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}49_60.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}61_72.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}73_84.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}85_96.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}97_102.grib2 >> ${dove}dati.csv

cut -d, -f1 --complement ${dove}dati.csv > ${dove}dati2.csv

cut -d, -f1 --complement ${dove}dati2.csv > ${dove}temp.csv

sed -i -e 's/val=//g' ${dove}temp.csv

go run ${dove}arpege.go > /var/www/html/luca/gfs/arpege.csv

---------------------------------------------------------------------------------

Confronto tra GFS NOAA e stazione meteo a terra

Confronto tra Arpege e stazione meteo a terra

Come si vede in valore assoluto i dati sono differenti ma e' comunque possibile creare una retta di correlazione tra i dati di previsione da terra ed i dati misurati dalla stazione meteo. Considerando che il primo dato e' un dato di previsione il coefficiente di correlazione non e' pessimo

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)

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