Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add rotated latitude longitude projection to mapping.py #3123

Merged
merged 5 commits into from
Jul 24, 2023
Merged

Add rotated latitude longitude projection to mapping.py #3123

merged 5 commits into from
Jul 24, 2023

Conversation

blaylockbk
Copy link
Contributor

@blaylockbk blaylockbk commented Jul 18, 2023

Description Of Changes

Checklist

The Canadian HRDPS model is on a rotated latitude/longitude projection grid. This PR adds the capability for MetPy to parse its projection.

Minimal working example with HRDPS data:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from pyproj import CRS
import pygrib
import cfgrib
import xarray as xr
import sys
import metpy

# First, you should download a HRDPS grib2 file from
# https://dd.weather.gc.ca/model_hrdps/continental/2.5km/00/00/
# I grabbed the TMP_AGL-2m file...
FILE = "20230722T00Z_MSC_HRDPS_TMP_AGL-2m_RLatLon0.0225_PT000H.grib2"
ds = xr.open_dataset(FILE)

# Parse the projparams with pygrib, because cfgrib doesn't do this for you
# (see https://github.com/ecmwf/cfgrib/issues/251)
with pygrib.open(str(FILE)) as grb:
    msg = grb.message(1)
    cf_params = CRS(msg.projparams).to_cf()

ds["gribfile_projection"] = None
ds["gribfile_projection"].attrs = cf_params
for var in list(ds):
    if var == "gribfile_projection":
        continue
    ds[var].attrs["grid_mapping"] = "gribfile_projection"

# Ok, now we are ready to check that the changes to metpy/mapping.py worked...
variables = [i for i in list(ds) if len(ds[i].dims) > 0]
ds = ds.metpy.parse_cf(varname=variables)
crs = ds.metpy_crs.item().to_cartopy()

# And plot the data...
ax = plt.subplot(111, projection=crs)
ax.coastlines()
ax.pcolormesh(ds.longitude, ds.latitude, ds.t2m, transform=ccrs.PlateCarree())
ax.set_global()
ax.gridlines()

image

@blaylockbk blaylockbk requested a review from a team as a code owner July 18, 2023 23:16
@blaylockbk blaylockbk requested review from dopplershift and removed request for a team July 18, 2023 23:16
@blaylockbk
Copy link
Contributor Author

@erin6541, I know you said you were going to work on this. Here is what I had so far, and it actually did work in my quick test.

@erin6541
Copy link

No worries, thanks for letting me know! Glad to hear your test worked.

@dopplershift
Copy link
Member

For testing, it's probably good to manually verify that the model data look right when plotted with the generated projection, but for adding to metpy just adding a more simple version in test_mapping.py, like what's there for other projections, should be sufficient.

@dopplershift dopplershift added Type: Enhancement Enhancement to existing functionality Area: Projections Pertains to projecting coordinates between coordinate systems labels Jul 24, 2023
@dopplershift dopplershift added this to the September 2023 milestone Jul 24, 2023
@dopplershift
Copy link
Member

Looks good! Thanks for the contribution!

@dopplershift dopplershift merged commit d49d6b3 into Unidata:main Jul 24, 2023
32 checks passed
@dopplershift
Copy link
Member

Congrats on your first merged PR to MetPy! 🎉 Hope to see you around in the future!

@blaylockbk
Copy link
Contributor Author

Thanks! I've been a MetPy user for years. Happy to have finally contributed something.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Projections Pertains to projecting coordinates between coordinate systems Type: Enhancement Enhancement to existing functionality
Projects
None yet
Development

Successfully merging this pull request may close these issues.

metpy/plots/mapping.py >> Unhandled projection: rotated_latitude_longitude
3 participants