Enhance the MET library code to read Rotated Lat/Lon data from CF-compliant NetCDF files.

[dwfncar]

NOTE: This has been identified as critical to Met Office deliverables - please charge 2799991

This request came from Steven Chan at the UK Met office via MET-Help:


   https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=87015

He sent us two files, a rotated lat/lon data file and one with variable lat/lon spacing near the poles. The goal is to make the former work in MET but not the latter since MET doesn't support variable mesh grids. These files can be found on dakota:
  /d3/projects/MET/MET_Help/chan_data_20180921/schan_eurocordex.nc

Here's the portion of the CF-convention document for Rotated Lat/Lon:
   http://cfconventions.org/cf-conventions/cf-conventions.html#_rotated_pole

It specifies the *NORTH POLE* instead of the south pole, as in GRIB2. Also, the grid dimension and spacing must be determined from the coordinate info.

[MET-1055] created by johnhg

[JohnHalleyGotway]

Discussed the high priority of this issue on 1/20/22 during the MetOffice developer's telecon. The MetOffice will provide additional NetCDF files for testing, and the MET team will ensure that it can be read.

[scole-mo]

All, please find the additional NetCDFs for testing below.

india_rotated_unrotated_grid.nc.gz

[mpm-meto]

Just a quick update from this end. From what I can gather now, when dealing with rotated grids for point observation applications it is best to compute the rotated coordinates of the observation before extracting the correct grid point(s) from the rotated model grid. This is because if you unrotate the grid first you no longer have a rectilinear grid and the whole thing just becomes a whole lot more awkward all around.

[JohnHalleyGotway]

The code on this feature branch is now able to read Rotated Lat/Lon Grids from CF-compliant NetCDF files.

For example, running plot_data_plane to plot it:

plot_data_plane \
india_rotated_unrotated_grid.nc \
orig.ps \
'name="number_of_lightning_flashes"; level="(*,*)";'

Results in:

Also running regrid_data_plane to regrid it:

regrid_data_plane india_rotated_unrotated_grid.nc \
"latlon 250 220 -10 54 0.25 0.25" regrid.nc \
-field 'name="number_of_lightning_flashes"; level="(*,*)";' \
-method MAX -width 5

And plotting the result yields:

Note that the features line up even though the shape of the underlying map data has changed.

Also note that I updated the logging so that at verbosity level 4, both the grid information and summary of the data plane values are printed:

DEBUG 4: Grid information:
DEBUG 4:    Projection: Rotated Lat/Lon
DEBUG 4:    Nx: 1200
DEBUG 4:    Ny: 1200
DEBUG 4:    rot_lat_ll: -24.000
DEBUG 4:    rot_lon_ll: 24.000
DEBUG 4:    delta_rot_lat: 0.040
DEBUG 4:    delta_rot_lon: 0.040
DEBUG 4:    wrapLon: false
DEBUG 4:    true_lat_south_pole: -72.000
DEBUG 4:    true_lon_south_pole: -85.000
DEBUG 4:    aux_rotation: 0.000
DEBUG 4: 
DEBUG 4: Data plane information:
DEBUG 4:    plane min: 0
DEBUG 4:    plane max: 18
DEBUG 4:    valid time: 20190321_010000
DEBUG 4:    lead time: 010000
DEBUG 4:    init time: 20190321_000000
DEBUG 4:    accum time: 010000

[JohnHalleyGotway]

@rgbullock I'm looking for a little advice. For this issue, we are enhancing the vx_data2d_nccf library code to parse Rotated Lat/Lon projections from CF-compliant NetCDF files.

This line of the feature_1055_read_rot_latlon branch, just hard-codes the "auxiliary rotation" value for the grid to 0. And I'm wondering whether or not we should be?

Here's a snapshot of the NetCDF header from a sample data file we have:

	int rotated_latitude_longitude ;
		rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ;
		rotated_latitude_longitude:longitude_of_prime_meridian = 0. ;
		rotated_latitude_longitude:earth_radius = 6371229. ;
		rotated_latitude_longitude:grid_north_pole_latitude = 72. ;
		rotated_latitude_longitude:grid_north_pole_longitude = 265. ;
		rotated_latitude_longitude:north_pole_grid_longitude = 0. ;

I'm wondering if the "aux_rotation" should be defined relative to the "longitude_of_prime_meridian" value? But when looking through the CF-conventions documentation, I only see examples where "longitude_of_prime_meridian = 0".

What do you recommend? Should I just leave the rotation hard-code to a value of 0? Or should I define it relative to the longitude_of_prime_meridian? If the latter, any advice on how?

[rgbullock]

Hi John - If I recall correctly, that auxiliary rotation parameter was just a placeholder in case we ever needed it. I think just having it default to zero would be fine. Randy

[JohnHalleyGotway]

@rgbullock thanks for confirming.

@jprestop please proceed with your review of PR #2041 as-is. I don't need to make any changes to it.

[mpm-meto]

@JohnHalleyGotway quick question on this... will this enable rotated grids in all tools? including GridDiag?

[JohnHalleyGotway]

@mpm-meto yes, the substantive changes in the corresponding pull request are to nccf_file.cc in the vx_data2d_nccf library. All MET tools call the same library code to read gridded data files. So yes, I fully expect it should work just find in GridDiag. I haven't explicitly tested that, but I don't anticipate any problems. Please let me know if you find otherwise!

[mpm-meto]

I've tested plot_data_plane on forecast and observed fields and it produced postscripts fine. I've also run GridDiag on these files and produced a netcdf file as output.