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)

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

NOAA 18

Una volta capito come funziona ci si prende gusto. A questo giro ho provato a ricevere i dati dal satellite NOOA18, un satellite meteo in orbita polare

Il passaggio sull'Europa e' avvenuto alle 21:15 (GMT+1) del 28 gennaio (fonte https://www.heavens-above.com/)

La caratteristica di questo satellite e' che trasmette sui 137.912 MHz con una ampiezza di banda di 38 KHz. In WebSDR per allargare l'ampiezza di banda si deve selezionare FM e poi con il mouse selezionare ciascun estremo della banda e trascinarlo fino a quando non viene coperta tutta l'estenzione del segnale. Il segnale e' caratterizzato da una serie di righe parallele nella modalita' waterfall

Waterfall

Spettro

Vista la frequenza di esercizio non potevo usare il WebSdr utilizzato per la  ISS e mi sono spostato su questo link http://erc-websdr.esa.int/

Di fatto il problema principale con NOOA e' possedere una antenna a polarizzazione circolare destra

Il file .wav deve essere processato attraverso il software WxtoImg ma e' posssibile utilizzare anche https://noaa-apt.mbernardi.com.ar/

Se si usa WxtoImg si deve prima prendere il file .wav in uscita da WebSDR e processarlo con Audacity per ricampionarlo da 8000 Hz s 11025 Hz. Si salva il file ma cio' comporta che il nuovo file abbia un timestamp diverso dal file iniziale. Per correggere questo problema si utilizza il comando

touch -a -m -t 202001282130.00 t1.wav

dove il timestamp e' formato YYYYMMDDHHMM.SS

Il software si chiama in questo modo perche' il protocollo di trasmissione e' denominato WeFax ed e' in pratica molto simile al buon vecchio sistema Fax che si usava in ufficio

Il servizio APT e' invece una trasmissione di due immagini relative a due bande (in questo caso nel termico) di 2080 x 909 pixels

Elaborazioni con WxToImg

Temperatura Mare

MCIR

Tenpeartura

se si usa noaa-apt i passi di cambio ora e resampling sono eseguiti in automatico e questo e' il risultato. Di fatto con questo software non possono essere eseguite elaborazioni. Le due bande vengono visualizzate a fianco con inclusa la banda di sincronia del segnale 

Per confronto questa e' la stessa immagine acquisita da https://www.osservatoriometeoesismicoperugia.it/satellitipolari/index.html

La differenza principale e' che qui le nuvole sono bianche mentre in wxtoimg risultano grigie