reading rotated latlon netcdf files

[perthsb]

Hi,
I am wondering is there any way MET/MET+ can read rotated latlon grid in netcdf file (like dynf*.nc or phyf*.nc from RRFS) model. I tried plot_data_plane, but it gives error like not a valid data file. Attached is ncdump of that netcdf file.
dump.txt
It can however read UPP processed rotated latlon grib2 files, attached is wgrib2 from a sample file
dump_grib2.txt

Thanks
Partha

[jprestop]

Hi @perthsb. Can you please post a log file or all of the output, including the command that you ran? Any additional information with regard to input files, machine name, etc. would be beneficial as well?

[perthsb]

@jprestop on Dogwood I tried following command
plot_data_plane /lfs/h2/emc/stmp/partha.bhattacharjee/Test39/stmp/2024021700/fcst_fv3lam/phyf000.nc test.ps 'name="ext550"; level="(0,64,,)";' -v 5
DEBUG 1: Opening data file: /lfs/h2/emc/stmp/partha.bhattacharjee/Test39/stmp/2024021700/fcst_fv3lam/phyf000.nc
DEBUG 4: Met2dDataFileFactory::new_met_2d_data_file() -> created new Met2dDataFile object of type "FileType_None".
ERROR :
ERROR : plot_data_plane -> file "/lfs/h2/emc/stmp/partha.bhattacharjee/Test39/stmp/2024021700/fcst_fv3lam/phyf000.nc" not a valid data file
ERROR :

plot_data_plane /lfs/h2/emc/ptmp/partha.bhattacharjee/Test39/ptmp/prod/rrfs.20240217/00/rrfs.t00z.prslev.f000.grib2 test2.ps 'name="AOTK"; level="A0";' -v 5 generates .ps file

[DanielAdriaansen]

Hi @perthsb. @jprestop was able to grab the phyf000.nc file you are using. In order to read NetCDF files with MET, they need to be one of three formats:

1. CF-compliant NetCDF
2. NetCDF output from another MET tool
3. NetCDF created by running the p_interp or wrf_interp utility on WRF output.

Looking at the requirements for CF-compliance in the MET user's guide here:
https://met.readthedocs.io/en/latest/Users_Guide/data_io.html#requirements-for-cf-compliant-netcdf

I was able to note right away some missing required info in your file, but I thought I would try and force MET to read it anyway. Here is the command I used:

plot_data_plane /lfs/h2/emc/stmp/partha.bhattacharjee/Test39/stmp/2024021700/fcst_fv3lam/phyf000.nc test.ps 'name="ext550"; level="(0,64,,)"; file_type=NETCDF_NCCF;' -v 5

And here is the output:

DEBUG 1: Opening data file: ./phyf000.nc
DEBUG 4: Met2dDataFileFactory::new_met_2d_data_file() -> created new Met2dDataFile object of type "FileType_NcCF".
DEBUG 4: NcCfFile::open() -> parsing units for the time variable "hours since 2024-02-17 00:00:00"
DEBUG 4: parse_cf_time_string() -> parsed NetCDF CF convention time unit string "hours since 2024-02-17 00:00:00"
DEBUG 4: 		 as a reference time of 20240217_000000 and 3600 second(s) per time step.
DEBUG 5: add_to_unixtime() -->20240217_000000 plus 0.01 times 3600 seconds = 20240217_000036
DEBUG 4: NcCfFile::open() -> could not extract init time from the "forecast_reference_time" variable.
DEBUG 4: NcCfFile::open() -> could not extract init time from file name.
DEBUG 4: NcCfFile::get_grid_from_dimensions() -> The coordinate variable "nchars" does not exist.
ERROR  : 
ERROR  : NcCfFile::read_netcdf_grid() -> Couldn't figure out projection from information in netCDF file.
ERROR  : 

Unfortunately, I don't think this file meets all of the requirements to use directly with MET. You can however, use Python embedding to read this file. Let us know if you'd like more help with that.

[perthsb]

Thanks @DanielAdriaansen for your answer. Can you provide help on how to read this with Python embedding.

[DanielAdriaansen]

Sure. Here is a script to read the file:

import xarray as xr

# Read netcdf file with Xarray
ds = xr.open_dataset('phyf000.nc')

# Which variable and level is desired
varname = 'ext550'
level = 64

# Obtain just the 2D variable 
# (select the level and remove any size 1 dimensions, e.g. time)
met_data = ds[varname].isel(pfull=level).squeeze()

# Dictionary of grid attributes
grid_attrs = {}
grid_attrs['type'] = "Rotated LatLon"
grid_attrs['name'] = "rrfs"
grid_attrs['rot_lat_ll'] = -37.0
grid_attrs['rot_lon_ll'] = 299.0
grid_attrs['delta_rot_lat'] = 0.025
grid_attrs['delta_rot_lon'] = 0.025
grid_attrs['Nlat'] = 2961
grid_attrs['Nlon'] = 4881
grid_attrs['true_lat_south_pole'] = -35.0
grid_attrs['true_lon_south_pole'] = 247.0
grid_attrs['aux_rotation'] = 0.0

# Dictionary of attributes required by MET
attrs = {}
attrs['valid'] = "20240217_000000"
attrs['init'] = "20240217_000000"
attrs['lead'] = "000000"
attrs['accum'] = "000000"
attrs['name'] = varname
attrs['long_name'] = 'extinction'
attrs['level'] = ""
attrs['units'] = "ext_units"
attrs['grid'] = grid_attrs

# Obtain a NumPy N-D array from the Xarray object and 
# reorder the data from S-N to N-S
met_data = met_data.values[::-1]

And the command to call it with plot_data_plane:

plot_data_plane PYTHON_NUMPY test.ps 'name="read_rrfs_netcdf.py";' -v 5`

And the resulting graphic that I get:

Does the graphic look reasonable to you? Let us know if you have any questions about the Python code!

[perthsb]

Thanks @DanielAdriaansen Graphic looks good and same as I get from Grib2 file (which is UPP processed but rotated latlon file) generated from the same netcdf file used here.

[hsoh-u]

set_attr_grid is for this case.
The grid definition of rotlatlon (from Appendix B Map Projections, Grids, and Polylines):

rotlatlon Nx Ny lat_ll lon_ll delta_lat delta_lon true_lat_sp true_lon_sp aux_rotation

Here is a command with set_attr_grid at seneca:

/d1/projects/MET/MET_regression/develop/NB20241108/MET-develop/bin/plot_data_plane  \
  /d1/projects/METplus/discussions/2507/phyf000.nc test.ps \
  'name="ext550"; level="(0,64,*,*)"; file_type=NETCDF_NCCF; set_attr_grid="rotlatlon 4881 2961 -37.0 299.0 0.025 0.025 -35.0 247.0 0.0";' \
  -v 4