DHOB (IU5SGN): NetCDF4

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