Hi,
I would like to share my experience on the visualization of ECMWF netcdf analysis with vapor. It works but it certainly could be improved so feel free to comment!
I convert both surfaces and 3D parameters.
First, my python script (see below) that convert ECMWF netcdf file in a new netcdf format more suitable for VAPOR conversion. Mainly, the latitudes are flipped, and the level as well. Not that an ELEVATION parameter seems to be required, I computed by other means and appended it to the ECMWF netcdf source file.
Then, the creation and population of the vdf file is easy through ncdfvdfcreate and ncdf2vdf (see below).
My only concern is about the handling of missing value / fill value.
As you can see in the python script, I use the fill_value argument in the createVariable method.
I tried many things, including using default_fillvals method of netCDF4, in combination with the -miss_attr option of ncdfvdfcreate but any values are interpreted by vaporgui as valid.
Eventually, I fixed the fill_value at 1.10 the maximum value and I play with the histogramm opacity in DVR tab of vaporgui, but it's not totally satisfactory.
#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
import netCDF4 as nc
from netCDF4 import Dataset,default_fillvals
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
#
# Read surfaces outputs from ECMWF
#-------------------------------------------------------------------------
fname = 'classes/filtre_SURF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc'
hdle = nc.Dataset(fname)
# lat/lon coord
lat_coord = hdle.variables['latitude'][:]
lon_coord = hdle.variables['longitude'][:]
time_coord = hdle.variables['time'][:]
# time,nlev,lat,lon
sp = hdle.variables['sp'][:,:,:]
st = hdle.variables['stl1'][:,:,:]
t2m = hdle.variables['t2m'][:,:,:]
d2m = hdle.variables['d2m'][:,:,:]
u10 = hdle.variables['u10'][:,:,:]
v10 = hdle.variables['v10'][:,:,:]
#
# Write new format
#-------------------------------------------------------------------------
fname_out = 'classes/VAPOR_filtre_SURF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc'
hdle_out = nc.Dataset(fname_out,'w')
hdle_out.createDimension("lat",lat_coord.size)
hdle_out.createDimension("lon",lon_coord.size)
hdle_out.createDimension("time",time_coord.size)
fval = default_fillvals['u8']
#hdle_out._FillValue = fval
#hdle_out.missingvalue = fval
lat = hdle_out.createVariable("lat",'f',("lat"))
lat.long_name="Latitude"
lat.units="degrees_north"
lat._CoordinateAxisType = "Lat";
# WARNING: latitudes are flipped
lat[:] = lat_coord[::-1]
lon = hdle_out.createVariable("lon",'f',("lon"))
lon.long_name="Longitude"
lon.units="degrees_east"
lon._CoordinateAxisType = "Lon";
lon[:] = lon_coord[:]
time = hdle_out.createVariable("time",'i',("time"))
time.long_name="time"
time.units="hours since 1900-01-01 00:00:0.0"
time._CoordinateAxisType = "Time";
time[:] = time_coord[:]
# WARNING: latitudes are flipped
fval = st.max()*1.1
stout = hdle_out.createVariable('stl1', 'f',("time","lat","lon"),fill_value=fval)
stout[:,:,:] = st[:,::-1,:]
stout.missing_value = fval
fval = t2m.max()*1.1
tout = hdle_out.createVariable('t2m', 'f',("time","lat","lon"),fill_value=fval)
tout[:,:,:] = t2m[:,::-1,:]
fval = d2m.max()*1.1
dout = hdle_out.createVariable('d2m', 'f',("time","lat","lon"),fill_value=fval)
dout[:,:,:] = d2m[:,::-1,:]
fval = u10.max()*1.1
uout = hdle_out.createVariable('u10', 'f',("time","lat","lon"),fill_value=fval)
uout[:,:,:] = u10[:,::-1,:]
fval = v10.max()*1.1
vout = hdle_out.createVariable('v10', 'f',("time","lat","lon"),fill_value=fval)
vout[:,:,:] = v10[:,::-1,:]
hdle_out.close()
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
#
# Read profiles outputs from ECMWF
#-------------------------------------------------------------------------
fname = 'classes/filtre_PROF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc'
hdle = nc.Dataset(fname)
# lat/lon coord
lat_coord = hdle.variables['latitude'][:]
lon_coord = hdle.variables['longitude'][:]
time_coord = hdle.variables['time'][:]
level_coord = hdle.variables['level'][:]
# time,nlev,lat,lon
t = hdle.variables['t'][:,:,:,:]
q = hdle.variables['q'][:,:,:,:]
clwc = hdle.variables['clwc'][:,:,:,:]
# REQUIRED but not a direct ouput from ecmwf
ELEV = hdle.variables['ELEVATION'][:,:,:,:]
#
# Write new format
#-------------------------------------------------------------------------
fname_out = 'classes/VAPOR_filtre_PROF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc'
hdle_out = nc.Dataset(fname_out,'w')
# dimensions
hdle_out.createDimension("lat",lat_coord.size)
hdle_out.createDimension("lon",lon_coord.size)
hdle_out.createDimension("time",time_coord.size)
fval = default_fillvals['u8']
#hdle_out._FillValue = fval
#hdle_out.missingvalue = fval
# Warning: latitudes are flipped
lat = hdle_out.createVariable("lat",'f',("lat"))
lat.long_name="Latitude"
lat.units="degrees_north"
lat._CoordinateAxisType = "Lat";
lat[:] = lat_coord[::-1]
lon = hdle_out.createVariable("lon",'f',("lon"))
lon.long_name="Longitude"
lon.units="degrees_east"
lon._CoordinateAxisType = "Lon";
lon[:] = lon_coord[:]
time = hdle_out.createVariable("time",'i',("time"))
time.long_name="time"
time.units="hours since 1900-01-01 00:00:0.0"
time._CoordinateAxisType = "Time";
time[:] = time_coord[:]
# sub sampling of levels
subsample = 1
levmax = 60
sublev = level_coord[:levmax:subsample]
hdle_out.createDimension("hybrid",sublev.size)
level = hdle_out.createVariable("hybrid",'f',("hybrid"))
level.long_name="Levels"
level._CoordinateAxisType = "GeoZ";
level._CoordinateZisPositive = "down";
level.units = 'sigma'
level[:] = sublev[:]
#print "Liste des niveaux"
#print level[:]
nlev = level.size
# Warning: latitudes are flipped
# Warning: levels are flipped
fval = t.max()*1.1
tout = hdle_out.createVariable('t', 'f',("time","hybrid","lat","lon"),fill_value=fval)
for klev,lev in enumerate(np.arange(nlev)[:levmax:subsample]):
tout[:,nlev-klev-1,:,:] = t[:,lev,::-1,:]
fval = q.max()*1.1
qout = hdle_out.createVariable('q', 'f',("time","hybrid","lat","lon"),fill_value=fval)
for klev,lev in enumerate(np.arange(nlev)[:levmax:subsample]):
qout[:,nlev-klev-1,:,:] = q[:,lev,::-1,:]
fval = clwc.max()*1.1
clout = hdle_out.createVariable('clwc', 'f',("time","hybrid","lat","lon"),fill_value=fval)
for klev,lev in enumerate(np.arange(nlev)[:levmax:subsample]):
clout[:,nlev-klev-1,:,:] = clwc[:,lev,::-1,:]
fval = ELEV.max()*1.1
elout = hdle_out.createVariable('ELEVATION','f',("time","hybrid","lat","lon"),fill_value=fval)
for klev,lev in enumerate(np.arange(nlev)[:levmax:subsample]):
elout[:,nlev-klev-1,:,:] = ELEV[:,lev,::-1,:]
hdle_out.close()
# NCDFCREATE
ncdfvdfcreate -vdc2 -missattr _Fillvalue -mapprojection "+proj=latlong +ellps=sphere" -extents -15400000:1100000:0:-11000000:4000000:16000 -timedims time classes/VAPOR_filtre_PROF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc classes/VAPOR_filtre_SURF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc ecmwf.vdf
# POPULATE
ncdf2vdf -timedims time -vars stl1:t2m:d2m:u10:v10 classes/VAPOR_filtre_SURF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc ecwmf.vdf
ncdf2vdf -timedims time -vars t:q:clwc classes/VAPOR_filtre_PROF01_0d25_20150601_00-06-12-18_region_Lon.-140.-100_Lat.10.40.nc ecmwf.vdf
Can you post an 'ncdump -h' of one of your source NetCDF files?