Transforming NARR LCC coordinate system to Lat Lon

[rubenjbehnke]

I have tried to find a way to transform the Lambert Conic Conformal projection (meters) in NARR data to WGS 84, with no success whatsoever. These files have only 'x' and 'y' as coordinates, so I can't get lat_lon_grid_deltas and other lat lon dependent functions to work. If anyone could help me out, that would be wonderful. Here is some code to get the data I am working with. Also, how do I extract data for a certain location, such as 44.325, -96.785 from this grid the right way? I can get a lot of things to work, but in a hacky sort of way, and am finding documentation for what seems like basic operations to be quite lacking:)

import numpy as np
import xarray as xr
import pandas as pd
from siphon.catalog import TDSCatalog
import fnmatch
import metpy.calc as calc
from metpy.units import units

yr = 1979
narr_data = TDSCatalog('https://rda.ucar.edu/thredds/catalog/files/g/ds608.0/3HRLY/' + str(yr) + '/catalog.xml')
ds = narr_data.datasets[np.min(np.where(fnmatch.filter(narr_data.datasets, '*3D*'))):np.max(np.where(fnmatch.filter(narr_data.datasets, '*3D*')))+1]

d = 0
lots_of_data = xr.open_dataset(ds[d].access_urls['OPENDAP']).metpy.parse_cf()

lots_of_data['u-component_of_wind_isobaric'].metpy.cartopy_crs

Out[132]: 
<Derived Projected CRS: +proj=lcc +ellps=sphere +a=6367470.0 +b=6367470.0  ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Lambert Conic Conformal (1SP)
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich

I am relatively new to Python, so I suppose that is adding to my struggles here. Otherwise, thank you!

[DanielAdriaansen]

Hi,

Thanks for sharing your code, it was helpful. You are close. This line:
lots_of_data = xr.open_dataset(ds[d].access_urls['OPENDAP']).metpy.parse_cf()

actually adds a Coordinate Reference System (CRS) to your lots_of_data variable, and this is the key to unlocking the rest of the puzzle. Here's what that looks like:

>>> lots_of_data.metpy.parse_cf().metpy_crs
<xarray.DataArray 'metpy_crs' ()>
array(<metpy.plots.mapping.CFProjection object at 0x7fc9ecdf9790>,
      dtype=object)
Coordinates:
    metpy_crs  object Projection: lambert_conformal_conic

From the documentation here:
https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#coordinates-and-coordinate-reference-systems

there is a part about your situation here:
https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#coordinates-and-coordinate-reference-systems:~:text=Finally%2C%20there%20are,the%20accessor%20documentation.

So, you can obtain the 2D latitude/longitude values for each grid cell in the NARR Lambert Conformal Conic grid like this:

data_with_lat_lon = xr.open_dataset(ds[d].access_urls['OPENDAP']).metpy.parse_cf().metpy.assign_latitude_longitude()
print(data_with_lat_lon.coords)
Coordinates:
  * x          (x) float32 -5.629e+06 -5.597e+06 ... 5.635e+06 5.668e+06
  * y          (y) float32 -4.61e+06 -4.577e+06 ... 4.317e+06 4.35e+06
  * time       (time) datetime64[ns] 1979-01-01 ... 1979-01-03T21:00:00
    reftime    (time) datetime64[ns] 1979-01-01 ... 1979-01-03T21:00:00
  * isobaric   (isobaric) float32 100.0 125.0 150.0 175.0 ... 950.0 975.0 1e+03
    metpy_crs  object Projection: lambert_conformal_conic
  * isobaric1  (isobaric1) float32 850.0
  * isobaric2  (isobaric2) float32 600.0 650.0 700.0 725.0 ... 950.0 975.0 1e+03
    latitude   (y, x) float64 1.0 1.104 1.208 1.312 ... 47.16 46.88 46.59 46.3
    longitude  (y, x) float64 -145.5 -145.3 -145.1 ... -2.633 -2.559 -2.486

I'm not familiar with WGS 84, but to properly locate the point 44.325, -96.785 in the NARR data, you ideally should compute the great circle distance (i.e. tunnel distance) between EVERY point in the grid and 44.325, -96.785, since it is not a latitude/longitude grid but rather a Lambert Conformal Conic projection. Taking the difference between all the latitudes and longitudes and 44.325, -96.785 will get you close (and perhaps be correct near the center of the grid), however it will not return accurate matches the further you get away from the center of the grid because of how the grid is draped over the Earth.

Here is resource I've used for understanding this problem that details several approaches. The NARR grid is not terribly high resolution so you may be able to get away with looping instead of taking the KDTree approach. The equations and concepts here though are still applicable to finding the distance:
https://www.unidata.ucar.edu/blogs/developer/en/entry/accessing_netcdf_data_by_coordinates

Hope this is helpful!

[rubenjbehnke]

Thank you! Works great.