Issue plotting WRF Data with Lambert Projection -Coastlines overlay not being displayed on top of data

Question

I recently started using python via jupyter on windows to plot netCDF data. However I am having issues plotting output data I have from a WRF run that uses Lambert Conformal projection. I cannot seem to get the coastlines to be overlayed and lined up with the data plot i.e. the coastlines overlay either does not show up at all or is slightly offset from the data.Could this be due to the Lambert projection of the data causing issues?

Here is the link to access the WRF output nc file I was using: https://drive.google.com/file/d/1s-0RLcaqIG0OoFEZtsE6_1xsE9B0DM1_/view?usp=sharing

The output of the code I used is attached.

I also don't seem to be able to use 'Dataset' on the WRF output files I am using, eg. ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00"). This normally results in the following error:

61 lons = wrfnc.variables[lonvar][:]
    663 # Need to check all times
--> 664 for i in py3range(lats.shape[-3]):
    665     start_idxs = [0] * len(lats.shape)  # PyNIO does not support ndim
    666     start_idxs[-3] = i

IndexError: tuple index out of range

I am quite new to all this so apologies if these are silly questions.

Below is the code I have tried to use. Any assistance in resolving this errors would be appreciated.

from netCDF4 import Dataset
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

ds = xr.open_dataset("output_mean.nc")
temp2 = ds.T2[0,:,:]
temp2

fig = plt.figure(figsize=(5, 5))
ax1 = plt.axes(projection=ccrs.LambertConformal())
ax1.coastlines()

temp2.plot(ax=ax1, cmap='jet', transform=ccrs.LambertConformal())

ax1.set_extent([-140, -75, 17, 45], crs=ccrs.PlateCarree())

ax1.coastlines()

plt.show()

Answer 1

There are a few different issues I think. It might also be more clear to avoid using the Xarray Dataset for the plotting. It's usually fine for convenience, but it also obfuscates some steps that might add to the confusion in this case. It also seems like Xarray is not correctly picking up on the coordinate dimensions, which would probably require some modifications if you want to use the Xarray Dataset for plotting directly.

I would start by reading the required data from the Xarray object to Numpy arrays:

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

with xr.open_dataset("output_mean.nc") as ds:
    ds = ds.isel(Time=0) # strip the time dimension
    lons = ds["XLONG"].to_numpy()
    lats = ds["XLAT"].to_numpy()
    data = ds["T2"].to_numpy()

The projection of the coordinates aren't specified as far as I can tell, but "WGS84" or ccrs.PlateCarree() in Cartopy, is often a good guess.

The default LambertConformal from Cartopy has a different definition than the data in your dataset, you can use the values defined in the dataset to specify the same projection, it looks like it should be something like this:

map_proj = ccrs.LambertConformal(
    central_latitude=39, 
    central_longitude=-101, 
    standard_parallels=(32, 46),
)

If you compare both you'll notice that the definition above doesn't "warp" the data, compared to the default definition. If the data is on a regular grid, plotting it in the same projection avoids having to do any interpolation when resampling to the different projection. It depends on the requirements whether that's needed or not.

One thing that you do incorrect is specifying LambertConformal as the transform of the data. But with Cartopy you should specify the projection of the source data (from which you want to transform), it should match the projection of the coordinates that you provide, which in this case are the lat/lon coordinate arrays. The target projection, where you want to transform to, is defined by the projection of the axes.

fig, ax = plt.subplots(
    figsize=(8, 5), facecolor="w",
    subplot_kw=dict(projection=map_proj ),
)

cm = ax.pcolormesh(lons, lats, data, cmap='turbo', transform=ccrs.PlateCarree())
cb = fig.colorbar(cm, ax=ax, shrink=.5)

# ax.set_extent([-140, -75, 17, 45], crs=ccrs.PlateCarree())
ax.coastlines()

If you omit setting the extent, it will show all data in the plot as the default:

Using the default definition of LambertConformal() makes it look like:

It's of course perfectly valid to plot in a different projection than the input data, but it's good to be aware when doing so.

An alternative would be to plot with ax.imshow (assuming it's a regular grid), that allows you to plot the data as it is (without resampling), but it would require knowing the extent of the data. It's probably possible to get that by reprojecting the corner lat/lon values to the projection as defined in the dataset.