Add support for the Lambert Azimuthal Equal Area grids

[michelleharrold]

Describe the New Feature

As of 10/25/22 - Identified as CRITICAL for Met Office:

MET assumes a spherical Earth. The EUMETNET (radar) dataset is on a Lambert Azimuthal Equal Area grid, which uses an ellipsoidal Earth. Randy already has some code for latlon_to_xy. The xy_to_latlon code still needs to be implemented, and the code needs to be implemented in MET (and tested accordingly).

Acceptance Testing

Same data is located on Cheyenne:
/gpfs/fs1/p/ral/jntp/AF/gswr/eumetnet (in HDF5 format)

This data can be moved to a more convenient location. Some sample data has already been shared with Randy.

Tests will need to be done to make sure the Lambert Azimuthal Equal Area projection is working as anticipated (will work with Randy to establish necessary tests).

Time Estimate

Estimate the amount of work required here.
Issues should represent approximately 1 to 3 days of work.

Sub-Issues

Consider breaking the new feature down into sub-issues.
No sub-issues required.

Relevant Deadlines

Goal is to target the met-10.1.0_beta1 release scheduled for 5/13/21.

Funding Source

Define the Metadata

Assignee

• Select engineer(s) or no engineer required: Randy B
• Select scientist(s) or no scientist required: Michelle H

Labels

• Select component(s)
• Select priority
• Select requestor(s)

Projects and Milestone

• Review projects and select relevant Repository and Organization ones or add "alert:NEED PROJECT ASSIGNMENT" label
• Select milestone to next major version milestone or "Future Versions"

Define Related Issue(s)

Consider the impact to the other METplus components.

• METplus, MET, METdatadb, METviewer, METexpress, METcalcpy, METplotpy
  No impacts.

New Feature Checklist

See the METplus Workflow for details.

• Complete the issue definition above, including the Time Estimate and Funding source.
• Fork this repository or create a branch of develop.
  Branch name: feature_<Issue Number>_<Description>
• Complete the development and test your changes.
• Add/update log messages for easier debugging.
• Add/update unit tests.
• Add/update documentation.
• Push local changes to GitHub.
• Submit a pull request to merge into develop.
  Pull request: feature <Issue Number> <Description>
• Define the pull request metadata, as permissions allow.
  Select: Reviewer(s), Project(s), Milestone, and Linked issues
• Iterate until the reviewer(s) accept and merge your changes.
• Delete your fork or branch.
• Close this issue.

[JohnHalleyGotway]

From @michelleharrold:

The deliverable for this project is due 31 May, so the sooner the better!

For reference - the EUMETNET data is what this is needed for. It is currently in HDF5, so the goal is to have you write the python code to convert to MET-friendly NetCDF and then have MET be able to handle the lambert azimuthal equal area grid. I suppose we could also have you include a regridding in the python script to bypass needing the lambert azimuthal equal area grid...but that seems like a cop out? But there is a hard deadline of 31 May -- and that is when the analysis needs to be done, so we need about 2 weeks to actually run MET and analyze output. Is this reasonable?

[JohnHalleyGotway]

From @rgbullock

I finished up the laea grid code. I've also tested it and it seems to be working well. I decided to go with the series expansion rather than the iterative scheme in the code, because it turns out that the iteration doesn't always converge.

The definition for the grid Michelle and Tracy are using is contained in the source file "laea_grid_defs.h". It's currently called tracy_data ... you'll probably want to change the name to something else.

For the benefit of those who will ultimately have the task of maintaining my grid code, the main reference for this work is

"Map Projections --- A Working Manual", by John P. Snyder
U. S. Geological Survey Professional Paper 1395

I'll let someone else write the code for wiring this up to the 2D data factory class. It should just be a matter of populating the LaeaData struct from the info in the data file.

I guess that's it unless there's something I've forgotten about.

One thing: Since this grid is on a non-spherical earth, any code that uses it will need to link to libvx_geodesy.a when compiling, so someone should probably add that to any of the Makefile.am files where it's needed so that autoconf knows about it.

[JohnHalleyGotway]

@rgbullock has finished development for this support. Added @JohnHalleyGotway to the assignee list so he can work on integrating these changes into the rest of the MET codebase.

[JohnHalleyGotway]

Here is a link to the GRIB2 definition of Lambert Azimuthal Equal Area grids:
https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-140.shtml

Along with an ECMWF grid template:
https://apps.ecmwf.int/codes/grib/format/grib2/templates/3/140

And here is a sample LAEA GRIB2 file I found on the MetOffice website:
https://metoffice.apiconnect.ibmcloud.com/metoffice/development/projection

wgrib2 -V ukv_agl_temperature_1.5_12.grib
	Lambert Azimuthal Equal Area grid: (1042 x 970) input WE:SN output WE:SN res 48
	Lat1 44.517153 Lon1 2164.600777 Cen Lon -2.500000 Std Par 54.900000
	Dx 2000.000000 m Dy 2000.000000 m mode 48

@rgbullock, how can we translate these GRIB2 projection definitions to the entries in the LaeaData struct that you've already added?

[JohnHalleyGotway]

@rgbullock can you please provide an update on this issue? Are you able to update the MET library code to instantiate a LAEA grid from a GRIB2 file?

[JohnHalleyGotway]

Thanks @robdarvell, I just grabbed those files but did not find what I expected:

ncdump -h 202211150600_nimrod_ng_radar_rainrate_composite_1km_UK.nc
	float rainfall_rate_composite(projection_y_coordinate, projection_x_coordinate) ;
...
		rainfall_rate_composite:grid_mapping = "transverse_mercator" ;
	int transverse_mercator ;
		transverse_mercator:grid_mapping_name = "transverse_mercator" ;
		transverse_mercator:longitude_of_prime_meridian = 0. ;
		transverse_mercator:semi_major_axis = 6377563.396 ;
		transverse_mercator:semi_minor_axis = 6356256.91 ;
		transverse_mercator:longitude_of_central_meridian = -2. ;
		transverse_mercator:latitude_of_projection_origin = 49. ;
		transverse_mercator:false_easting = 400000. ;
		transverse_mercator:false_northing = -100000. ;
		transverse_mercator:scale_factor_at_central_meridian = 0.9996012717 ;

Looks like this is for a Transverse Mercator projection rather than Lambert Azimuthal Equal Area.

Do you have any sample data for the latter that I can use to develop/test this issue?

[robdarvell]

Apologies @JohnHalleyGotway for the duff data. I was told the radar files were on the LAEA grid, but I hadn't actually checked through the ncdump before I sent them.

I have placed a new temperature and visibility file onto the ftp server, and checked that the ncdump talked about LAEA on them. Hopefully (!) they are what you need.
They cover the UKV (ie UK) domain.

Rob

[TaraJensen]

Trying to get this going again. @JohnHalleyGotway - did we ever grab the files @robdarvell put up on the FTP server? If not - hopefully @robdarvell can post them again.

[robdarvell]

I have put some alternate LAEA grid files on the ftp server in darvell_data
(they are probs of thresholds being met fields)