crear netcdf

# Cafe de programacion CIMA/DCAO/IFAECI

# Ejemplo de como crear un fichero netCDF con python3

# Librerias
import numpy as np
from netCDF4 import Dataset as NetCDFFile

# Parametros
##

# nombre fichero
ncfilen = 'test.nc'

# Dimensiones
dimx = 100
dimy = 50
dimz = 20
dimt = 8

# Datos de las proyeccion de los datos
ilon = -80.
elon = -40.
ilat = -30.
elat = -50.

# presiones
presvals = [1000., 975., 950., 925., 900., 850., 800., 750., 700., 650., 600., 550., \
  500., 450., 400., 350., 300., 250., 200., 100.]

# tiempo
## reftime
reftime = '1949-12-01 00:00:00'
timeu = 'hours'
## diferencia en timeu del 18 de Agosto de 2022, respecto el refdate de las unidades 
## de tiempo
diffbasetime = 637392

#######    #######
## MAIN
    #######

onc = NetCDFFile(ncfilen, 'w')

# Creando las dimensiones (siguiendo convenciones CF, https://cfconventions.org/)
newdim = onc.createDimension('lon', dimx)
newdim = onc.createDimension('lat', dimy)
newdim = onc.createDimension('pres', dimz)
newdim = onc.createDimension('height2', 1)
newdim = onc.createDimension('height10', 1)
newdim = onc.createDimension('time', None)

# Creando las variables dimensionales

# Asumimos proyeccion regular lon/lat (campos 1D de longitudes y latitudes)
newvar = onc.createVariable('lon', 'f', ('lon'))
dlon = (elon - ilon)/(dimx-1)
lons = np.arange(ilon, elon + dlon, dlon)
newvar[:] = lons
newvar.setncattr('standard_name', 'longitude')
newvar.setncattr('long_name', 'Longitude')
newvar.setncattr('units', 'degrees_east')

newvar = onc.createVariable('lat', 'f', ('lat'))
dlat = (elat - ilat)/(dimy-1)
lats = np.arange(ilat, elat + dlat, dlat)
newvar[:] = lats
newvar.setncattr('standard_name', 'latitude')
newvar.setncattr('long_name', 'Latitude')
newvar.setncattr('units', 'degrees_north')

newvar = onc.createVariable('height2', 'f', ('height2'))
newvar[:] = 2.
newvar.setncattr('standard_name', 'height')
newvar.setncattr('lon_name', 'Height above ground')
newvar.setncattr('units', 'm')
newvar.setncattr('positive', 'up')

newvar = onc.createVariable('height10', 'f', ('height10'))
newvar[:] = 10.
newvar.setncattr('standard_name', 'height')
newvar.setncattr('lon_name', 'Height above ground')
newvar.setncattr('units', 'm')
newvar.setncattr('positive', 'up')

newvar = onc.createVariable('pres', 'f', ('pres'))
newvar[:] = np.array(presvals)*100.
newvar.setncattr('standard_name', 'pressure')
newvar.setncattr('lon_name', 'Pressure')
newvar.setncattr('units', 'Pa')
newvar.setncattr('positive', 'down')

newvar = onc.createVariable('time', 'f', ('time'))
newvar[:] = np.arange(0, 24, 3) + diffbasetime
newvar.setncattr('standard_name', 'time')
newvar.setncattr('long_name', 'Time')
newvar.setncattr('units', timeu + ' since ' + reftime)
newvar.setncattr('calendar', 'gregorian')

# Variables

newvar = onc.createVariable('tas', 'f', ('time', 'height2', 'lat', 'lon'),           \
  fill_value=1.e20)

tas = np.full((dimt, 1, dimy, dimx), 305.)
newvar[:] = tas
newvar.setncattr('standard_name', 'air_temperature')
newvar.setncattr('long_name', 'temperature of the air')
newvar.setncattr('units', 'K')
newvar.setncattr('coordinates', 'time height2 lat lon')

newvar = onc.createVariable('uas', 'f', ('time', 'height10', 'lat', 'lon'),          \
  fill_value=1.e20)
uas = np.full((dimt, 1, dimy, dimx), 5.)
newvar[:] = uas
newvar.setncattr('standard_name', 'eastward_wind')
newvar.setncattr('long_name', 'eastward wind speed of the air')
newvar.setncattr('units', 'ms-1')
newvar.setncattr('coordinates', 'time height10 lat lon')

newvar = onc.createVariable('vas', 'f', ('time', 'height10', 'lat', 'lon'),          \
  fill_value=1.e20)
vas = np.full((dimt, 1, dimy, dimx), 3.)
newvar[:] = vas
newvar.setncattr('standard_name', 'northward_wind')
newvar.setncattr('long_name', 'northward wind speed of the air')
newvar.setncattr('units', 'ms-1')
newvar.setncattr('coordinates', 'time height10 lat lon')

newvar = onc.createVariable('zg', 'f', ('time', 'pres', 'lat', 'lon'),               \
  fill_value=1.e20)
zg = np.full((dimt, dimz, dimy, dimx), 0.)
# Para el ejemplo, el primer nivel esta por encima de los 1000 hPa
zg[:,0,:,0:dimx//2] = 1.e20
newvar[:] = zg
newvar.setncattr('standard_name', 'geopotential_height')
newvar.setncattr('long_name', 'geopotential height')
newvar.setncattr('units', 'm2s-2')
newvar.setncattr('coordinates', 'time pres lat lon')

# Atributos globales
onc.setncattr('author', 'anonymous')
onc.setncattr('institution1', 'Centro de ' + 'Investigaciones del Mar y la Atm' +    \
  chr(243) + 'sfera (CIMA)')
onc.setncattr('university', 'Universidad de Buenos Aires (UBA)')
onc.setncattr('institution2', 'Instituto Franco-Argentino sobre Estudios de Clima '+ \
  'y sus Impactos (CNRS, IRL-3351 IFAECI)')
onc.setncattr('version', '0.0')