WRF in Python Feature Comparison

[jthielen]

Here is the first of my promised write-ups from #2! It is a little messy, but hopefully is sufficient to be a good reference. I plan on having the initial "user stories" write-up ready before next week's meeting.

---

Motivation: How do features in existing purpose-built tools (wrf-python, salem) compare with those currently in more general purpose packages (xgcm, MetPy, etc.)? What are the gaps in the existing ecosystem in meeting user needs?

Comparison of general-purpose tools against wrf-python feature stack

The following goes through all wrf-python's features and compares them to what is available in MetPy (or other general-purpose packages). The goal is to identify what features are needed either in xwrf or added to MetPy/elsewhere to be able to supersede all current wrf-python features.

Diagnostics

wrf-python key	Description	Status Relative to MetPy
avo	Absolute Vorticity	Implemented: mpcalc.absolute_vorticity
eth/theta_e	Equivalent Potential Temperature	Implemented: mpcalc.equivalent_potential_temperature
cape_2d	2D CAPE (MCAPE/MCIN/LCL/LFC)	Implemented for single-profile only: mpcalc.cape_cin
cape_3d	3D CAPE and CIN	Implemented for single-profile only: mpcalc.cape_cin
ctt	Cloud Top Temperature	No equivalent (though, unclear how to make a useful generalized formulation)
cloudfrac	Cloud Fraction	No equivalent (though, unclear how to make a useful generalized formulation)
dbz	Reflectivity	Use directly from output
mdbz	Maximum Reflectivity	Use .max('bottom_top') on reflectivity DataArray
geopt/geopotential	Geopotential for the Mass Grid	Field needing WRF-specific processing
geopt_stag	Geopotential for the Vertically Staggered Grid	Field needing WRF-specific processing
helicity	Storm Relative Helicity	Implemented for single-profile only: mpcalc.storm_relative_helicity
lat	Latitude	Use directly from output
lon	Longitude	Use directly from output
omg/omega	Omega	Implemented: mpcalc.vertical_velocity_pressure
p/pres	"Full Model Pressure (in specified units)"	Field needing WRF-specific processing
pressure	Full Model Pressure (hPa)	Field needing WRF-specific processing
pvo	Potential Vorticity	Implemented: mpcalc.potential_vorticity_baroclinic
pw	Precipitable Water	Implemented for single-profile only: mpcalc.precipitable_water
rh	Relative Humidity	Implemented: mpcalc.relative_humidity_from_mixing_ratio
rh2	2m Relative Humidity	Implemented: mpcalc.relative_humidity_from_mixing_ratio
slp	Sea Level Pressure	Open PR Pending: Unidata/MetPy#1700
T2	2m Temperature	Field needing WRF-specific processing
ter	Model Terrain Height	Field needing WRF-specific processing
td2	2m Dew Point Temperature	Implemented, but multiple functions possibly needed depending on exact source fields (hopefully will be made easier with Unidata/MetPy#2070)
td	Dew Point Temperature	Implemented, but multiple functions possibly needed depending on exact source fields (hopefully will be made easier with Unidata/MetPy#2070)
tc	Temperature in Celsius	Use .metpy.convert_units('degC') or .pint.to('degC') on temperature DataArray
th/theta	Potential Temperature	Implemented: mpcalc.potential_temperature
temp	Temperature (in specified units)	Field needing WRF-specific processing
tk	Temperature in Kelvin	Field needing WRF-specific processing
times	Times in the File or Sequence	Field needing WRF-specific processing
xtimes	"XTIME Coordinate (if applicable)"	Field needing WRF-specific processing
tv	Virtual Temperature	Implemented: mpcalc.virtual_temperature
twb	Wet Bulb Temperature	Implemented for single-profile only: mpcalc.wet_bulb_temperature
updraft_helicity	Updraft Helicity	Open Issue: Unidata/MetPy#2213
ua	U-component of Wind on Mass Points	Field needing WRF-specific processing
va	V-component of Wind on Mass Points	Field needing WRF-specific processing
wa	W-component of Wind on Mass Points	Field needing WRF-specific processing
uvmet10	"10 m U and V Components of Wind Rotated to Earth Coordinates"	Field needing WRF-specific processing (also, PyProj/cartopy vector transforms could use some helpers with xarray)
uvmet	"U and V Components of Wind Rotated to Earth Coordinates"	Field needing WRF-specific processing (also, PyProj/cartopy vector transforms could use some helpers with xarray)
wspd_wdir	"Wind Speed and Direction (wind_from_direction) in Grid Coordinates"	Implemented: mpcalc.wind_speed and mpcalc.wind_direction
wspd_wdir10	"10m Wind Speed and Direction (wind_from_direction) in Grid Coordinates"	Implemented: mpcalc.wind_speed and mpcalc.wind_direction
uvmet_wspd_wdir	"Wind Speed and Direction (wind_from_direction) Rotated to Earth Coordinates"	Partially implemented: mpcalc.wind_speed and mpcalc.wind_direction (PyProj/cartopy vector transforms could use some helpers with xarray)
uvmet10_wspd_wdir	"10m Wind Speed and Direction (wind_from_direction) Rotated to Earth Coordinates"	Partially implemented: mpcalc.wind_speed and mpcalc.wind_direction (PyProj/cartopy vector transforms could use some helpers with xarray)
z/height	Model Height for Mass Grid	Field needing WRF-specific processing
height_agl	Model Height for Mass Grid (AGL)	Field needing WRF-specific processing
zstag	Model Height for Vertically Staggered Grid	Field needing WRF-specific processing

Interpolation

wrf-python routine	Status Relative to MetPy
wrf.interplevel()	Doable, but too complecated...should be a job for xgcm
wrf.vertcross()	mpinterp.cross_section or mpinterp.interpolate_to_slice
wrf.interpline()	mpinterp.cross_section or mpinterp.interpolate_to_slice
wrf.vinterp()	Doable, but too complecated...should be a job for xgcm

Lat/Lon <-> XY Routines

These packages have different approaches to coordinate handling. Ideally, xwrf should provide both grid projection space xy dimension coordinates and auxillary lat lon coordinates, with the WRF-specific coordinate metadata translated to CF conventions so that MetPy, xgcm, and cf-xarray can interpret it and nothing special needs to be handled in-package.

Mapping Helper Routines

Should be obviated by good dataset metadata interpretable by xarray.

Meeting User Needs

Here's the list of the "core features" needed/requested in a Pangeo-centric WRF package:

• Pythonic, xarray-centric API
• Full Pangeo stack integration, particularly Dask for parallelization and lazy-loading
• Large collection of performant diagnostic routine calculations that remain portable (i.e., use numba rather than cython or f2py if vectorized NumPy code is insufficient) and are automatically available (i.e., Dask delayed)
• Support for both automatic destaggering (for typical analysis workflows) and working on staggered grid (with xgcm)
• CF compliant metadata with geospatial dimension coordinates (i.e., Salem-like, as needed for integration with other packages like MetPy, cf-xarray, xesmf, etc.)
• Community-oriented development philosophy (e.g.; integrate with existing packages where functionality exists rather than rolling our own; any new, non-WRF-specific functionality needed here gets contributed upstream then reused here)

Where do these features exist currently and where are further developments needed?

xarray-centric API

Currently

• wrf-python falls short on this front due to getvar()-routine-based, and reliance on netcdf-python datasets, rather than xr.Dataset-based API
  • Also, poorly compatible metadata
• Vanilla xarray also falls short, as WRF's coordinates are a bit too esoteric to be properly handled natively
• salem has one option for a good API, but is limited in scope and pending deprecation
• xwrf has a promising approach, which is a good target to build around

Unmet Needs

Compatibility layer to

• translate WRF output metadata/coordinates to CF (for all variables present in output that are not otherwise handled below)
  • this includes west_east, south_north, and bottom_top dimension coordinates
• munge "raw" split/offset fields to meaningful physical variables (core variables like pressure, geopotential, potential temperature)
• optional (but perhaps by default) automatic destaggering of staggered variables
• do all this in a lazy-loaded/Dask-leveraged way
• present as a xr.Dataset through xarray backend and/or xarray accessor API

These are the essential goals for xwrf!

Pangeo stack integration (particularly Dask)

Currently

• Dask support often only through one-off/custom side tools

Unmet Needs

• See above goals for xwrf

Diagnostic routines

Currently

• wrf-python has its list, but poorly integrates with other packages like MetPy if they are insufficient (e.g., frontogensis, divergence), and perhaps has too many for a dataset-specific tool
• MetPy has large suite, as is fitting for a general-purpose toolkit, but is substantially hampered on large datasets for anything parcel-based
• salem has very limited variables

Upcoming Developments

• MetPy is being funded to improve performance of calculations, particularly parcel-based ones so that they are grid-capable

Unmet Needs

• A tool that can map hopefully all, but at least most, variables in WRF output to either
  • CF metadata for the variable as-is
  • Munged with other fields to provide a physically meaningful variable (that in turn can be described well with CF)
• xarray-friendly geospatial vector transform utilities

Staggered and Unstaggered Grids

Currently

• Salem and wrf-python have good support for destaggering variables
• Poor metadata for working with WRF output on its true (staggered) grid

xwrf goals

• Follow lead of salem and wrf-python and provide easy destaggered variable access (though, do it with Dask arrays automatically)
• Map the metadata to CF/COMODO conventions so that xgcm can handle all staggered grid needs

CF compliant metadata with geospatial dimension coordinates

Covered above...but see salem for example of what to provide, and MetPy for what is needed

[kmpaul]

@jthielen @lpilz @michaelavs @dopplershift @LejoFlores @jukent @andersy005 @erogluorhan @mgrover1 @bsu-wrudisill:

Please chime in on this discussion! When the plans get solidified, perhaps we can move this to a design document (GitHub Project?) that we can all modify and share.

[jukent]

This is a great outline of the work ahead of us. Do I understand correctly that all the wrf-python keys with the note "Field needing WRF-specific processing" are priority?

[lpilz]

As I understood it, this just means that we have to build some infrastructure for these keys ourselves rather than being able to rely on metpy infrastructure or contribute to it.

[singledoggy]

Is it possible to add a temporary middleware to place the computation results of wrf-python into the xwrf data?"

[mgrover1]

Not sure what is meant here?

[singledoggy]

It seems that xwrf can eventually replace wrf-python, but many diagnostic variables in wrf-python are still useful in certain situations. Currently, when these diagnostic variables are needed, we still have to extract them using wrf-python's getvar() method and manually merge them into xwrf's data. Is there a simpler way to achieve a more consistent reading method without explicitly using getvar(), while still using getvar() in the background?

For me, using getvar() is convenient, but for beginners, learning multiple libraries adds a lot of unnecessary hassle. This hinders me from recommending xwrf to them, even though I believe xwrf's integration with dask is a significant advantage.

[mgrover1]

This would be a pretty significant level of effort here - especially since much of the getvar() functionality is running fortran/ncl routines under the hood - are you interested in making a contribution here?

[singledoggy]

Thank you for pointing that out. I tried starting with the basic temperature variables, but it seems that in wrf-python, this also calls Fortran's dcomputetk, and some of the cloud diagnostics are more complex. I think fully migrating the calculation of diagnostic variables to pure Python is very labor-intensive, and incorporating Fortran calculations might completely change the computational process. I'm beginning to understand why there hasn't been much progress on this work; it might also be beyond my capabilities. However, is it possible for WRF to read the diagnostic variables calculated by wrf-python and incorporate the coordinate system in the form of xwrf? I'm not quite sure if doing so makes sense.

What I mean is, if there are already some netCDF files or variables calculated by wrf-python, based on this premise, we could just change the data structure. Perhaps we could further add delayed computation with something like dask.delayed. I think that wrf-python can also be wrapped by Dask, although it might not be very efficient.

I apologize for my confusing explanation. In most cases, I am just a user and haven't designed a large project from scratch.