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:
|NCL||not that high||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-(xlim_d01-xlim_d01)/15, xlim_d01+(xlim_d01-xlim_d01)/15]) ax.set_ylim([ylim_d01-(ylim_d01-ylim_d01)/15, ylim_d01+(ylim_d01-ylim_d01)/15]) # d01 box ax.add_patch(mpl.patches.Rectangle((xlim_d01, ylim_d01), xlim_d01-xlim_d01, ylim_d01-ylim_d01, fill=None, lw=3, edgecolor='blue', zorder=10)) ax.text(xlim_d01+(xlim_d01-xlim_d01)*0.05, ylim_d01+(ylim_d01-ylim_d01)*0.9, 'D01', size=15, color='blue', zorder=10) # d02 box ax.add_patch(mpl.patches.Rectangle((xlim_d02, ylim_d02), xlim_d02-xlim_d02, ylim_d02-ylim_d02, fill=None, lw=3, edgecolor='black', zorder=10)) ax.text(xlim_d02+(xlim_d02-xlim_d02)*0.05, ylim_d02+(ylim_d02-ylim_d02)*1.1, 'D02', size=15, color='black', zorder=10) # d03 box ax.add_patch(mpl.patches.Rectangle((xlim_d03, ylim_d03), xlim_d03-xlim_d03, ylim_d03-ylim_d03, fill=None, lw=3, edgecolor='red', zorder=10)) ax.text(xlim_d03+(xlim_d03-xlim_d03)*0.1, ylim_d03+(ylim_d03-ylim_d03)*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.