Visualizing WRF domain

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
NCL not that high namelist.wps
wrf-python high WRF output
mine 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[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:

  1. My domain uses Lamber Conformal projection (lcc), which causes extra trouble in showing the lat/lon gridlines in the figure.
  2. 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.).
  3. 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.
  4. 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.

10 thoughts on “Visualizing WRF domain

  1. Hi,

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

    Thanks,
    stella

    • Hi,

      You might want to preprocess the ETOPO1 data a little bit. The size of raw global ETOPO1 is too huge, so I regridded it to 0.1 degrees, and only clipped the North America region.

      • Hi, I downloaded the ETOOPO1 data,and could you please tell me how to regrid it too 0.1 and convert specified range to nc format ?

      • Hi,

        I used the CDO remapbil to regrid the nc format data to 0.1 degree. Regarding the format conversion, I was using gdal_translate, and the command was “gdal_translate -of nc “. Hope this helps.

  2. Hi,

    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!
    Xin

    • Hi Xin,

      I think the issue is from that the ETOPO1 file you downloaded is not in the current directory. Do a “ls ETOPO1_Ice_g_gdal.grb” and see if you can get it. You may want to use the full path to the file if it is not in the current directory, something like “/home/xxx/Download/ETOPO1_Ice_g_gdal_grb”.

  3. Hello,
    How could I make ax1.add_feature(OCEAN, edgecolor=’k’, facecolor=’lightblue’) work? Since I am using cmap=”terrain” ocean color is first color in this cmap.

    Thanks,
    Yandy

    • Hi Yandy,

      Sorry I did not quite get your question. Are you trying to use the first color from “terrain” cmap as the OCEAN color? If that is the case, you can find the specific color code of the first color following this discussion (https://stackoverflow.com/questions/25408393/getting-individual-colors-from-a-color-map-in-matplotlib).

      If you mean the “lightblue” is overwritten by your “terrain” cmap, then you just need to manually set the layer orders using “zorder”. It goes from 0 to say 10, and you just give the layer you want to put on top a larger number. Say you give the “terrain” cmap layer zorder=2, and OCEAN layer as zorder=3 (i.e., ax1.add_feature(OCEAN, edgecolor=’k’, facecolor=’lightblue’, zorder=3)). This way OCEAN is always above your “terrain” layer. Let me know if this solves your question.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s