There are many ways to visualize the spatial domain of WRF. I have used several of them in my previous papers, but they turned out to be not so elegant for publications. Therefore, while writing the dissertation recently, I developed some functions that can directly digest WPS namelist (rather than WRF output in some cases) to derive the domain boundaries. Figure 1 is an example plot.

Figure 1. Visualization of WRF domain from home-made Python. Background is the topography from ETOPO1 dataset.

Before this one, I had also tried NCL and wrf-python packages. See below for their outputs:

Figure 2. Visualization of WRF domain from NCL (plotgrids.ncl).

Figure 3. Visualization of WRF domain from wrf-python.

Here is a quick comparison of these 3 approaches:

Method Customization Input not that high namelist.wps wrf-python WRF output namelist.wps

It is obvious that NCL can provide a quick overview of the domain configuration, which is super useful in setting up the model (when iterative modifications to the domain are needed). But it does require significant modification to the script to make the plot a better fit for publications. Also, by default, it presents the D01 domain as the whole figure, which is sometimes a little confusing. I am sure we can modify this, but I did not put too much effort here.

wrf-python, on the other hand, provides easy to use functions that easily get the projection information, domain boundaries from WRF output. Below are the codes for Figure 3.

import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import netCDF4 as nc from cartopy import crs from cartopy.feature import NaturalEarthFeature import wrf def get_plot_element(infile): rootgroup = nc.Dataset(infile, 'r') p = wrf.getvar(rootgroup, 'RAINNC') #lats, lons = wrf.latlon_coords(p) cart_proj = wrf.get_cartopy(p) xlim = wrf.cartopy_xlim(p) ylim = wrf.cartopy_ylim(p) rootgroup.close() return cart_proj, xlim, ylim infile_d01 = 'domain_plot/wrfout_d01.2010-05-01_00:00:00.nc' cart_proj, xlim_d01, ylim_d01 = get_plot_element(infile_d01) infile_d02 = 'domain_plot/wrfout_d02.2010-05-01_00:00:00.nc' _, xlim_d02, ylim_d02 = get_plot_element(infile_d02) infile_d03 = 'domain_plot/wrfout_d03.2010-05-01_00:00:00.nc' _, xlim_d03, ylim_d03 = get_plot_element(infile_d03) fig = plt.figure(figsize=(10,8)) ax = plt.axes(projection=cart_proj) states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', name='admin_1_states_provinces_shp') ax.add_feature(states, linewidth=0.5) ax.coastlines('50m', linewidth=0.8) # d01 ax.set_xlim([xlim_d01[0]-(xlim_d01[1]-xlim_d01[0])/15, xlim_d01[1]+(xlim_d01[1]-xlim_d01[0])/15]) ax.set_ylim([ylim_d01[0]-(ylim_d01[1]-ylim_d01[0])/15, ylim_d01[1]+(ylim_d01[1]-ylim_d01[0])/15]) # d01 box ax.add_patch(mpl.patches.Rectangle((xlim_d01[0], ylim_d01[0]), xlim_d01[1]-xlim_d01[0], ylim_d01[1]-ylim_d01[0], fill=None, lw=3, edgecolor='blue', zorder=10)) ax.text(xlim_d01[0]+(xlim_d01[1]-xlim_d01[0])*0.05, ylim_d01[0]+(ylim_d01[1]-ylim_d01[0])*0.9, 'D01', size=15, color='blue', zorder=10) # d02 box ax.add_patch(mpl.patches.Rectangle((xlim_d02[0], ylim_d02[0]), xlim_d02[1]-xlim_d02[0], ylim_d02[1]-ylim_d02[0], fill=None, lw=3, edgecolor='black', zorder=10)) ax.text(xlim_d02[0]+(xlim_d02[1]-xlim_d02[0])*0.05, ylim_d02[0]+(ylim_d02[1]-ylim_d02[0])*1.1, 'D02', size=15, color='black', zorder=10) # d03 box ax.add_patch(mpl.patches.Rectangle((xlim_d03[0], ylim_d03[0]), xlim_d03[1]-xlim_d03[0], ylim_d03[1]-ylim_d03[0], fill=None, lw=3, edgecolor='red', zorder=10)) ax.text(xlim_d03[0]+(xlim_d03[1]-xlim_d03[0])*0.1, ylim_d03[0]+(ylim_d03[1]-ylim_d03[0])*0.8, 'D03', size=15, color='red', zorder=10) ax.set_title('WRF nested domain setup (2010Nash event)', size=20) plt.show() fig.savefig('WRF_domain_pywrf.big.png', dpi=600)

One can further simplify this code, as the major use of wrf-python here is just to create the projection (as in cart_proj = wrf.get_cartopy(p) ). The rest functions can all be implemented using just netCDF4 (or xarray) with some projection transformation, i.e. cart_proj.transform_point(….). Using Python also provides greater flexibility to improve the plot with other packages/data. One trouble with wrf-python, as can be seen above, is that it can only get information of one domain from a single file. Thus it requires several inputs when you have nested domains. This requies the WRF simulation (at least for the very first timestep) to be done before we can visualize the domain.

That’s the motivation of my home-made Python. Basically I would like to be able extract projection/boundary information from namelist.wps (just like NCL), but also own greater control on the final plot (so Python is preferred). I put the code I used for Figure 1 on Github, and here are some explanations:

My domain uses Lamber Conformal projection (lcc), which causes extra trouble in showing the lat/lon gridlines in the figure.

For now it only includes lcc projection. But you can easily find the logic of creating and switching to the other projections allowed in WPS (Mercator, pole, lat/lon, etc.).

Once you get the basic map to work, you can overlay whatever data (land use, topography, precipitation, …) you would like to show. This would make it more than just a domain configuration and a better fit for papers/slides.

For now, I only consider the situation where the ref_lat and ref_lon of D01 domain is the center of the D01 domain in the projected map. This is often the case of WPS setup, but it is possible (and it is allowed) that ref_lat and ref_lon are away from center. I may add this logic later.

I am putting the codes on Github ( https://github.com/lucas-uw/WRF-tools/blob/master/WRF_input_tools/Visualize_WPS_domain.ipynb ), and you are free to take it and modify it. In the meanwhile, if I have time I will polish it further and make it more user-friendly.

I downloaded the ETOPO1 data from this website:

https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/

Then you can convert the file to netCDF format or whatever you are comfortable with. If you shoot me an email (you can find it on my “official” website), I can send the converted netCDF file to you.

Reply

I try to use your codes to draw figure 1, but it shows “memory error”… How many memory you used?

Thanks,
stella

Reply

How did you convert grb to netcdf?
I’ve tried ‘cdo -f nc copy ETOPO1_Ice_g_gdal.grb ETOPO1_Ice_g_gdal.nc’, but I got this error:

cdo copy: Open failed on >ETOPO1_Ice_g_gdal.grb<
No such file or directory

Could you help me with this?

Thank you!
Reply

Privacy & Cookies: This site uses cookies. By continuing to use this website, you agree to their use.
To find out more, including how to control cookies, see here: Cookie Policy