lunedì 22 agosto 2022

Evento Stromboli 12 agosto su ERA5

 Il 12 agosto 2022 e' avvenuto un inteso evento meteo a Stromboli ..volevo vedere se era possibile ricavare i dati di precipitazione anche dal servizio ERA5 di Copernicus per confrontarli con la verita' a terra data dal pluviometro INGV

per scaricare i dati ho utilizzato il seguente script Python (dati del solo 12 agosto per il solo parametro Total Precipitation)

//////////////////////////////////

import cdsapi

c = cdsapi.Client()

c.retrieve(
'reanalysis-era5-single-levels',
{
'product_type': 'reanalysis',
'variable': 'total_precipitation',
'year': '2022',
'month': '08',
'day': '12',
'time': [
'00:00', '01:00', '02:00',
'03:00', '04:00', '05:00',
'06:00', '07:00', '08:00',
'09:00', '10:00', '11:00',
'12:00', '13:00', '14:00',
'15:00', '16:00', '17:00',
'18:00', '19:00', '20:00',
'21:00', '22:00', '23:00',
],
'area': [
39, 15.25, 38.75,
15.5,
],
'format': 'netcdf',
},
'12.nc')

//////////////////////////////////

i dati sono stati acquisiti in formato NetCDF4 perche' il formato GRIB non sembra compatibile con Wgrib2 (forse si tratta di grib1)..per leggere i dati ho prima usato Qgis per verificare di avere centrato la finestra geografica di ricerca



e poi ho utilizzato il seguente script Python per estrarre i dati (non ho trovato niente di gia' pronto per una conversione da NetCDF4 a testo)

import netCDF4

precip_nc_file = '12.nc'
nc = netCDF4.Dataset(precip_nc_file, mode='r')

nc.variables.keys()

lat = nc.variables['latitude'][:]
lon = nc.variables['longitude'][:]
time_var = nc.variables['time']
dtime = netCDF4.num2date(time_var[:],time_var.units)
precip = nc.variables['tp'][:]

print("lat "+str(len(lat)))
print("lon "+str(len(lon)))
print("prec "+str(len(precip)))
print("time "+str(len(dtime)))

# coordinate Stromboli 38.79184084355282, 15.2152918018263


for s in range(len(dtime)):
for l in range(len(lat)):
for g in range(len(lon)):
print(dtime[s], end=''),
print(";", end='')
print(lat[l], end=''),
print(";", end='')
print(lon[g], end=''),
print(";", end='')
# il primo indice e' quello dell'ora
# il secondo indice e' la latitudine
# il terzo indice e' la longitudine
print(precip[s][l][g], end='\n\r')

questo e' il grafico generato con i dati di ERA5
questo per confronto sono i dati dal pluviometro di Stromboli gestito da INGV


(la scala dei tempi e' differente ed anche l'ora risulta differente a causa dei differenti fusi)

In conclusione il dato da satellite sottostima l'evento ma la forma generale del grafico e' raffrontabile



Debian e Ideapad 310 MIIX310

Mi sono comprato il tablet/convertibile Lenovo Ideapad 310 MIIX non perche' ne avessi bisogno ma per il costo estremamente ridotto (50 euro) e .... mi sono messo nei guai

Con Windows 10 la macchina e' completamente inutilizzabile ma me lo aspettavo....cio' che non aspettavo era di dover lottare con i driver per Linux

Una sintesi

con Ubuntu la scheda di rete Realtek8723BS funziona, non funziona audio e webcam (la webcam e' un modello molto particolare non USB ne' PCI ..e' una OV2680..comoda su Arduino per niente su PC per mancanza di supporto)




con Manjaro in versione live funziona praticamente tutto (tranne webcam)...una volta installato e' sparita la scheda di rete

Debian nemmeno a parlarne...il supporto alla Realtek 8723BS esisteva nel kernel fino a qualche anno fa ma e' stato rimosso. Ho provato ad installare Devuan per avere i driver gia' installati ...funziona l'audio ma niente wifi e webcam

alla fine pero' Debian e' molto rapida su una macchina dalle risorse cosi' limitate e mi sono adattato con una scheda Wifi USB

Un aspetto fastidioso e' che essendo un tablet il sistema operativo parte in modalita' portrait, per ruotare lo schermo e il touchscreen si possono utilizzare i comani

xrandr -o right   
xinput set-prop 'FTSC1000:00 2808:1015' 'Coordinate Transformation Matrix' 0 1 0 -1 0 1 0 0 1

 

 

 sdfsfd

sabato 13 agosto 2022

Distanza Realsense

Misura la distanza al centro dell'immagine

Per rendere la misura piu' stabile e' stato utilizzato un ring buffer


## License: Apache 2.0. See LICENSE file in root directory.
## Copyright(c) 2015-2017 Intel Corporation. All Rights Reserved.

###############################################
## Open CV and Numpy integration ##
###############################################

import pyrealsense2 as rs
import numpy as np
from numpy_ringbuffer import RingBuffer
import cv2


RingBufferSize = 50
r = RingBuffer(capacity=RingBufferSize, dtype=np.float16)

# Configure depth and color streams
pipeline = rs.pipeline()
config = rs.config()

# Get device product line for setting a supporting resolution
pipeline_wrapper = rs.pipeline_wrapper(pipeline)
pipeline_profile = config.resolve(pipeline_wrapper)
device = pipeline_profile.get_device()
device_product_line = str(device.get_info(rs.camera_info.product_line))

found_rgb = False
for s in device.sensors:
if s.get_info(rs.camera_info.name) == 'RGB Camera':
found_rgb = True
break
if not found_rgb:
print("The demo requires Depth camera with Color sensor")
exit(0)

config.enable_stream(rs.stream.depth, 640, 480, rs.format.z16, 30)

if device_product_line == 'L500':
config.enable_stream(rs.stream.color, 960, 540, rs.format.bgr8, 30)
else:
config.enable_stream(rs.stream.color, 640, 480, rs.format.bgr8, 30)

# Start streaming
pipeline.start(config)

try:
while True:

# Wait for a coherent pair of frames: depth and color
frames = pipeline.wait_for_frames()
depth_frame = frames.get_depth_frame()
color_frame = frames.get_color_frame()
if not depth_frame or not color_frame:
continue

# Convert images to numpy arrays
depth_image = np.asanyarray(depth_frame.get_data())
color_image = np.asanyarray(color_frame.get_data())
#print(depth_image[320,240])
r.append(depth_image[320,240])
print(np.mean(r))

# Apply colormap on depth image (image must be converted to 8-bit per pixel first)
depth_colormap = cv2.applyColorMap(cv2.convertScaleAbs(depth_image, alpha=0.03), cv2.COLORMAP_JET)

depth_colormap_dim = depth_colormap.shape
color_colormap_dim = color_image.shape

# If depth and color resolutions are different, resize color image to match depth image for display
if depth_colormap_dim != color_colormap_dim:
resized_color_image = cv2.resize(color_image, dsize=(depth_colormap_dim[1], depth_colormap_dim[0]), interpolation=cv2.INTER_AREA)
images = np.hstack((resized_color_image, depth_colormap))
else:
images = np.hstack((color_image, depth_colormap))

# Show images
cv2.namedWindow('RealSense', cv2.WINDOW_AUTOSIZE)
startpoint=(310,240)
endpoint=(330,240)
color=(255,0,0)
thickness=3
images = cv2.line(images,startpoint,endpoint,color,thickness)
startpoint=(320,230)
endpoint=(320,250)


images = cv2.line(images,startpoint,endpoint,color,thickness)
cv2.imshow('RealSense', images)
cv2.waitKey(1)

finally:

# Stop streaming
pipeline.stop()

mercoledì 3 agosto 2022

Archivi di dati meteo

 Il sito https://cds.climate.copernicus.eu/cdsapp#!/search?type=dataset permette di scaricare in modo gratuito dati storici a carattere ambientale come il dataset ERA5 (meteo dati dal 1959 all'attuale)


I dati si possono scaricare dall'interfaccia web ma si possono anche utilizzare le API Python

Per fare cio' ci si deve iscrivere al sito per ricevere la API KEY (da ricavare dal proprio user profile) e la si deve inserire nel file HOME/.cdsapirc 

url: https://cds.climate.copernicus.eu/api/v2
key: 123456:123445-abcd-efgh-12345
a questo punto si puo' installare la libreria

per creare il programma si puo' andare sulla pagina web del dataset, selezionare i dati che si intendono scaricare e poi selezionare in fondo alla pagina il pulsante Show API Request. Basta copiare ed incollare nel file programma

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

import cdsapi
c = cdsapi.Client()
c.retrieve(
    'reanalysis-era5-single-levels',
    {
        'product_type': 'reanalysis',
        'format': 'grib',
        'year': '2022',
        'month': '07',
        'day': '01',
        'time': [
            '00:00', '01:00', '02:00',
            '03:00', '04:00', '05:00',
            '06:00', '07:00', '08:00',
            '09:00', '10:00', '11:00',
            '12:00', '13:00', '14:00',
            '15:00', '16:00', '17:00',
            '18:00', '19:00', '20:00',
            '21:00', '22:00', '23:00',
        ],
        'variable': '2m_temperature',
    },
    'download.grib')

mercoledì 27 luglio 2022

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


lunedì 25 luglio 2022

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

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

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




Macchie solari

 Riprese con Celestron 70/700 con filtro in mylar...avevo comprato il filtro solare l'anno scorso quando l'attivita' solare era al minimo e non erano presenti macchie😊😊




lunedì 14 marzo 2022

Platformio IDF MPU6050

Ho voluto provare a fare un esperimento programmando un ESP32 con il framework IDF al posto di quello Arduino (esperimento che rimarra' solitario perche' alla fine non ho trovato grandi vantaggi)

Ho avuto diverse difficolta' ad usare il monitor seriale di Platformio per cui ho usato Minicom


Il codice (si tratta sostanzialmente dell'esempio della libreria con impostate lo porte SDA SCL a GPIO 21 e 22 nel file i2c.c alla riga 23 e 24) 



sabato 12 marzo 2022

Anatomia di un conducimetro Neomet 29D

In ufficio e' morto il conducimetro da cantiere. Una volta rimosso dall'inventario e' arrivato il momento di sezionarlo




All'interno si puo' osservare la presenza di un microcontrollore ATMega 103L (8 bit, Flash 128 Kb, EEPROM 4 Kb, SRAM 4Kb), un integrato MAX232CE (RS232 porta seriale) ed un display da 4.3 pollici 94V0. Il display e' saldato alla base 



Al di sotto si trova un ICL7135CPI, un convertitore analogico digitale a 4 digit 1/2 ed un DC/DC converter tc7662cpa



Manuale 


 



venerdì 11 marzo 2022

Protocol Buffers in Golang

 Protocol Buffers (protobuf) e' una libreria di Google per la serializzazione di dati. Con il termine serializzazione si intende il processo di tradurre una "data structure" (array, struttura, albero binario.....) in un file che conservi oltre che i dati anche la struttura. Una forma tipica di serializzazione sono i file in formato JSON

Per utilizzare protobuf con GoLang si parte dalla definizione della struttura in file .proto. Un semplice modello puo' essere il seguente

syntax = "proto3";

package main;
import "google/protobuf/timestamp.proto";


message Person {
string name = 1;
int32 id = 2;
string email = 3;
google.protobuf.Timestamp last_updated = 5;
}

da notare che il file .proto deve avere lo stesso package del codice che lo utilizzera'

i numeri di fianco ai vari campi sono univoci

il file .proto deve essere compilato tramite il compilatore protoc. Su Debian l'installazione del compilatore avviene mediante

apt install protobuf-compiler
apt install golang-goprotobuf-dev

si crea il file .proto nella stessa directory del sorgente Go e si compila con

protoc --go_out=. *.proto

questo genera un nuovo file con estensione .go. Se si apre il file con editor di testo si osserva che oltre alla struttura dei dati sono state create delle funzioni tipo helper che facilitano l'immissione e gestione della struttura dati

per gestire i dati all'interno del progetto Go, nel caso si usino tipo dati di uso comune tipo datetime gia' formalizzati in protobuf si deve scaricare la loro definizione per esempio

go get github.com/golang/protobuf/ptypes/timestamp

i dati di un protobuf possono essere salvati su un file ma questo sara' binario e non immediatamente human readable

package main

import (
"fmt"
"io/ioutil"
"log"

proto "github.com/golang/protobuf/proto"
)

// apt install protobuf-compiler

// per compilare il file proto
// protoc --go_out=. *.proto

// go get github.com/golang/protobuf/ptypes/timestamp

func main() {
dati_persona := &Person{
Name: "Luca",
Id: 51, //per poco ancora
}

data, err := proto.Marshal(dati_persona)
if err != nil {
log.Fatal("marshaling error: ", err)
}

fmt.Println(data)

out, err := proto.Marshal(dati_persona)
if err != nil {
log.Fatalf("Serialization error: %s", err.Error())
}
if err := ioutil.WriteFile("dati.bin", out, 0644); err != nil {
log.Fatalf("Write File Error: %s ", err.Error())
}
fmt.Println("Write Success")

//Read from file
in, err := ioutil.ReadFile("dati.bin")
if err != nil {
log.Fatalf("Read File Error: %s ", err.Error())
}
dati_persona2 := &Person{}
err2 := proto.Unmarshal(in, dati_persona2)
if err2 != nil {
log.Fatalf("DeSerialization error: %s", err.Error())
}

fmt.Println("Read Success")
fmt.Printf("Nome %s\n", dati_persona2.GetName())
}




ESP32 IDF con PlatformIo

 


Utilizzando il file platformio.ini si possono configurare alcune impostazioni del progetto come le librerie da importare (anche direttamente da GitHub), la porta seriale dell'Esp32 (in caso il riconoscimento automatico non funzioni)

impostazione seriale
upload_port = /dev/ttyUSB3

importazioni librerie (attenzione che il software non risolve le dipendenze, queste devono essere gestite  a mano)
lib_deps = https://github.com/gabrielbvicari/esp32-mpu6050
https://github.com/gabrielbvicari/esp32-i2c_rw


#include "mpu6050/mpu6050.h"


pio lib install


lunedì 7 marzo 2022

Github Desktop su Debian

Trovo utile usare su Windows Github Desktop (piu' che altro perche' gestisce bene  i fork e le pull request). Lo stesso software puo' essere installato anche su Debian  

wget -qO - https://packagecloud.io/shiftkey/desktop/gpgkey | sudo apt-key add -

sh -c 'echo "deb [arch=amd64] https://packagecloud.io/shiftkey/desktop/any/ any main" > /etc/apt/sources.list.d/packagecloud-shiftky-desktop.list'

apt-get update

apt install github-desktop



venerdì 4 marzo 2022

Recupero Sticky Notes

Un collega mi telefona dicendo di avere perso dati importanti salvati nei PostIt virtuali di Windows 10 con una applicazione chiamata Sticky Notes....francamente non sapevo nemmeno che Windows avesse una applicazione di sistema di questo tipo ma quando mi dice che il file dove i dati vengono salvati e' chiamato plum.sqlite gli dico che provo a dare un'occhiata (ho lavorato su qualche progetto usando Sqlite)


In particolare la nota incriminata non risulta essere stata cancellata ma a schermo presenta una serie di ++++++++. Aprendo il file con effettivamante nella tabella note si vedono solo 3 righe di cui una con soli segno +. Il controllo di integrita' del database e' passato senza problemi quindi non ci sono corruzioni di indici

Stavo rinunciando quando mi sono ricordato che i database di solito non cancellano i dati direttamente ma solo dopo un esplicito comando ...ho provato a fare un export in csv da DB Browser for Sqlite per vedere se oltre ai dati mostrati venivano esportate anche le righe cancellate...bingo.. il file di esportazione presentava i dati cancellati in formato umanamente leggibile




lunedì 28 febbraio 2022

Chrome Os Flex

 Ho provato ChromeOs Flex... diciamo una mezza delusione (pur sapendo che si tratta di una Unstable)


Si tratta sostanzialmente di CloudReady rimarchiato (in molti punti e' presente ancora la scritta CloudReady)

Per avviare il container Linux ho dovuto attivare tutti i flag Crostini in http://flags ed in ogni caso il sottosistema Linux e' molto lento, sono riuscito ad avviare solo le Linux Apps e non un Window Manager

ChromeOS puro funziona ma si vede nettamente che non avendo un hardware ottimizzato come sui ChromeBooks non e' un prestante

L'idea di poter continuare hardware obsoleto con ChromeOs Flex al momento attuale non e' ottimale..tanto vale usare una distro Linux con poche pretese hardware 


Misteri

Un collega mi telefona dicendomi che il suo HD esterno non funziona piu' sottolineando che sente dei rumori (...ma i rumori non vengono dalle casse del PC😀). Gli dico che puo' essere la testina del HD che sbatte a fine corsa 

Non proprio convinto il mio collega apre il case dell'HD e mi manda un video (screenshot nella foto superiore)...la testina effettivamente sbatte a fine corsa ma mi sembra che manchi qualcosa...tipo la testina dell'HD...gli chiedo se quando ha aperto il case e' caduto fuori qualcosa ma si dice sicuro che non c'era niente di rotto

Gli chiedo il modello e gli mando una foto di archivio da Google Photo in cui ovviamente si vede la testina superiore presente



  Il mistero della sparizione della testina dell'HD 





venerdì 25 febbraio 2022

Quadtree in Go

metto questo abbozzo di codice perche', nonostante non funzionante, puo' essere interessante

Ai tempi di Fractint su Dos c'e un algoritmo che divideva il piano di quadrati, se i punti estremi del quadrato avevano lo stesso valore allora tutto il quadrato assumeva lo stesso colore dell'insieme di Mandelbrot, in caso contrario il quadrato veniva a sua volta diviso in 4 quadrati e si reiterava

In questo caso ho cercato di fare lo stesso usando le goroutine  ed la ricorsione


package main

import (
    "fmt"
    "math"
    "math/cmplx"
    "reflect"
    "strconv"
    "sync"

    "github.com/fogleman/gg"
)

var wg sync.WaitGroup

// ATTENZION    E il nome delle variabili deve essere in uppercase
// per poter esportare la variabile e
type punto struct {
    Min_x float64
    Max_x float64
    Min_y float64
    Max_y float64
    Itera int64
}

var punti []punto // slice di struct
var transi punto
var mu sync.Mutex

func calcola(a float64, b float64) int {
    var k int
    c := complex(a, b)
    k = 0

    x := complex(0.0, 0.0)
    for k = 0; k < 25; k++ {
        x_a := (x * x)
        x_a = x_a + c
        x = x_a
        if cmplx.Abs(x) > 2 {
            return k
        }
    }
    return k
}

func dividi(min_x float64, max_x float64, min_y float64, max_y float64) {
    defer wg.Done()

    x1 := calcola(min_x, min_y)
    x2 := calcola(max_x, max_y)

    if (x1 != x2) && (max_x-min_x > 0.005) {

        mezzo_x := (max_x - min_x) / 2
        mezzo_y := (max_y - min_y) / 2
        wg.Add(4)
        go dividi(min_x, min_x+mezzo_x, min_y, min_y+mezzo_y)
        go dividi(min_x+mezzo_x, max_x, min_y, min_y+mezzo_y)
        go dividi(min_x, min_x+mezzo_x, min_y+mezzo_y, max_y)
        go dividi(min_x+mezzo_x, max_x, min_y+mezzo_y, max_y)

    } else {
        fmt.Println(strconv.FormatFloat(min_x, 'f', 6, 32) + ";" + strconv.FormatFloat(min_y, 'f', 6, 32) + ";" + strconv.Itoa(x1))

        s := reflect.ValueOf(&transi).Elem()
        s.Field(0).SetFloat(min_x)
        s.Field(1).SetFloat(max_x)
        s.Field(2).SetFloat(min_y)
        s.Field(3).SetFloat(max_y)
        s.Field(4).SetInt(int64(x1))

        mu.Lock()
        punti = append(punti, transi)
        mu.Unlock()
    }

}

func grafica() {
    dc := gg.NewContext(400, 400)

    dx := 3.0 / 400
    dy := 2.0 / 400

    for _, value := range punti {
        fmt.Printf("%d\n", value)
        min_x := math.Round(value.Min_x*dx) + 200
        max_x := math.Round(value.Max_x*dx) + 200
        min_y := math.Round((value.Min_y * dy)) + 200
        max_y := math.Round((value.Max_y * dy)) + 200
        dc.DrawRectangle(min_x, max_x, min_y, max_y)
        dc.SetRGB255(value.Itera, 0, 0)
        dc.Fill()

    }

    dc.SavePNG("mandelbrot.png")
}

func main() {
    wg.Add(1)
    go dividi(-2.0, 1.0, -1.0, 1.0)
    wg.Wait()
    grafica()
}

Grafica base in Go

 package main


import "github.com/fogleman/gg"

func main() {
    dc := gg.NewContext(1000, 1000)
    //dc.Push()
    dc.DrawCircle(500, 500, 400)
    dc.SetRGB255(0, 0, 0)
    dc.Fill()
   
    dc.SetHexColor("#008080")
    dc.DrawRectangle(10,10,40,30)
    dc.SetRGB255(255, 0, 0)

    dc.Fill()
   
    dc.SavePNG("out.png")
}

LLama3 Anita

A seguito di questo post ho provato a vedere ho provato a vedere cosa accadeva ad utilizzare un modello specifico per la lingua italiana in...