Plotting WRF data using python
March 16, 2021: I have since then wrote a new story with an updated script here Plotting WRF data using python (wrf-python and cartopy edition)
Recently I have discovered the unholy matrimony of two of my favorite things: the Weather Research Forecasting (WRF) model and Quantum GIS (QGIS).
WRF is a numerical Weather Prediction (NWP) that allows for the simulation of the atmosphere using globally available data. QGIS on the other hand is a Geographic Information Systems platform that allows for manipulation or display of geospatial data.
Last 9 November 2018, the GIS4WRF plugin for QGIS was released. It allows for the complete steps for WRF simulation in a more visual manner. With GIS4WRF you can do everything from downloading input data, pre-processing input data, running simulations, visualizing, and finally post-processing results.
Back then when I used to simulate for my thesis I had to do things through the command line where I had to use terminal magic to make things happen. GIS4WRF would allow for easy and streamlined WRF simulations through QGIS.
As a quick test to know whether I still know how to run my weather simulations, I ran a simulation of Super Typhoon Yolanda using GIS4WRF.
Now I discovered a road block. Post-processing WRF data using GIS4WRF wasn’t as easy as you would think it would be. I initially wanted to plot the wind map of Super Typhoon Yolanda complete with the wind barbs, but I ended up with this.
Or it’s also possible that I didn’t possess the necessary skills just yet, but as GIS4WRF is a relatively new plugin, it kept crashing on me.
So python to the rescue. Particularly the following packages: gdal, basemap, matplotlib, and numpy.
The basic flow is the following
- Read WRF netcdf output
- Calculate wind speed
- Plot
From those three simple steps we have this:
Some caveats.
As much as GIS4WRF tries their best to make things user-friendly, WRF is still an advanced model. The plugin in itself won’t offer default physics and dynamics for WRF, as they say that “this is meant to be science so you should know the options to choose” assuming that you know how to simulate using WRF to begin with.
Next, no models are perfect. Outputs churned out by models such as WRF have their own limitations and must still be subject to expert analysis. Whenever you see models like the above, please always try to understand the context.
Finally, for the code that you may plot WRF output for yourself:
from osgeo import gdal
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy
import os# Specify the base path
PATH_base = 'D:/gis4wrf/projects/yolanda/'
# Specify the folder where *.nc files can be found
PATH_nc = os.path.join(PATH_base, 'run_wrf')
# Specify desired output folder
PATH_out = os.path.join(PATH_base, 'out')# Read a single file
# Remember, WRF output are nc files
# (even if the extension is absent)
PATH_file = os.path.join(PATH_nc, 'wrfout_d01_2013-11-01_00_00_00')# Constants
# Read a single time slice
ts = 56
# barb size
bs = 3
windmax = 30
windmin = 0# Read Data
ds_lon = gdal.Open('NETCDF:"' + PATH_file + '":XLONG')
ds_lat = gdal.Open('NETCDF:"' + PATH_file + '":XLAT')
ds_u10 = gdal.Open('NETCDF:"' + PATH_file + '":U10')
ds_v10 = gdal.Open('NETCDF:"' + PATH_file + '":V10')# for lon and lat we can use any array. in this case we use the
# index -1
lon = ds_lon.ReadAsArray()[-1]
lat = ds_lat.ReadAsArray()[-1]
u10 = ds_u10.ReadAsArray()[ts]
v10 = ds_v10.ReadAsArray()[ts]# for drawing meridians
meridians = numpy.arange(110,160,10.)
parallels = numpy.arange(0, 25, 10.)
fig = plt.figure(figsize=(30, 25))ax = fig.add_subplot(223)
ax.set_title("Yolanda surface wind.")
m = Basemap(llcrnrlon=110,llcrnrlat=0,urcrnrlon=160,urcrnrlat=25, resolution='l')# For wind barbs
x, y = m(lon, lat)
yy = numpy.arange(0, y.shape[0],bs)
xx = numpy.arange(0, x.shape[1],bs)
points = numpy.meshgrid(yy, xx)# Plot
m.pcolor(x, y, numpy.sqrt(u10*u10 + v10*v10), vmin = windmin, vmax = windmax, cmap = cm.viridis)
m.barbs(x[points], y[points], u10[points], v10[points], barbcolor='#CCCCCC')
m.drawmeridians(meridians, labels=[1,0,0,1])
m.drawparallels(parallels, labels=[1,0,0,1])
m.drawcoastlines()
m.colorbar()# Save
fig.savefig(os.path.join(PATH_out, 'yolanda_' + str(ts).zfill(2) + '.png'), bbox_inches = 'tight')
The code above was written with the help of this blog post: http://geoexamples.blogspot.com/2013/09/reading-wrf-netcdf-files-with-gdal.html
Happy coding, folks. :)
