venerdì 22 luglio 2022

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

giovedì 24 marzo 2022

Stima distanze con Aruco Tags ed OpenCV

Ho voluto provare questo progetto

https://github.com/GSNCodes/ArUCo-Markers-Pose-Estimation-Generation-Python

ho provato ad usare la webcam del portatile per provare a stimare le distanze mediante i tag Aruco (sono concettualmente simili a QrCode ma piu' robusti per il riconoscimento a scapito del numero di dati inseriti...di fatto solo un codice ID)..altri progetti usano gli AprilTag

Attenzione : e' necessario Python3

Per generare i tag si usa il codice sottostante

python3 generate_aruco_tags.py --id 24 --type DICT_4X4_50 --output tags/ -s 230

id = codice inserito nel tag
type = formato del tag (ne esistono molti...si puo' vedere l'elenco nel file utils.py)
output = e' il folder dove sono salvati i file png dei tag
s = dimensioni del tag

Esempio di tag

Il passo successivo e' quello di crea i file di calibrazione della camera. Si usa una scacchiera (checkerboard)...ho fatto alcune prove ed alcune immagini non funzionano ...ho usato quella proposta da OpenCv in formato 9x6
Si devono fare fotografie in cui si vede la scacchiera posizionata a diversa distanza e con diverse angolazioni e si mettono in una sottodirectory 

Checkerboard

Una delle immagini di calibrazione

Dopo aver raccolto una ventina di immagini di calibrazione (nel mio caso salvato nel folder pose) si lancia

python3 calibration.py --dir './pose/' -s 0.023  -w 9 -t 6

se tutto va a buon fine saranno create i file della matrice di calibrazione e dei coefficienti di distorsione


A questo punto si puo' procedere con la stima di distanza

python3 pose_estimation.py --K_Matrix calibration_matrix.npy --D_Coeff distortion_coefficients.npy --type DICT_4X4_50



il codice e' leggermente differente dal progetto originario perche' di fatto in origine vengono calcolate le matrici rvec e tvec ma non viene effettuato il calcolo della distanza. Le modifiche possono essere ritrovate al mio GitHub https://github.com/c1p81/Aruco_distance




Debugger integrato ESP32S3

Aggiornamento In realta' il Jtag USB funziona anche sui moduli cinesi Il problema risiede  nell'ID USB della porta Jtag. Nel modulo...