Repository navigation
Generalize vorticity (and related functions) to allow for map projections #893
Description
Some MetPy functions, such as gradient and advection, allow the user to supply arrays of deltas rather than a single value, but vorticity (and its relatives) does not. Vorticity is more complicated since derivatives of the "map factors" are also required, but it should be doable.
Here is the necessary computation according to the GEMPAK documentation (Appendix B2):
Since the map factors (m_x and m_y) are equivalent to the "nominal" grid spacing (\partial x and \partial y in the above equation) divided by the actual delta, an additional "nominal_delta" argument to vorticity would allow for the map factors to be computed. Since the WRF model provides the map factors in its output, allowing the map factors (and the nominal grid spacing) as optional arguments to vorticity in lieu of deltas would be useful as well.
Activity
dopplershift commentedon Jul 25, 2018
So while the docs for
vorticitycurrently say that dx/dy have to befloat(oops, that needs fixing), with the new finite differencing in 0.7, we can now handle vorticity for arbitrarily-spaced points. If you've tried to pass an array in for dx/dy, and it failed, that's a bug.As far as the mapping/projection issues are concerned, that's still an open issue (covered in #174). So I'm a little unclear--is the map factor there to handle irregularly-spaced points or to handle issues from projection? or both?
sgdecker commentedon Jul 27, 2018
OK, I see that
vorticitydoes compute something when I provide arrays for dx/dy.Regarding the map factors, they are fundamentally there to handle issues from projection. One of those issues is that the data array is irregularly spaced, but another issue is that computing vorticity (more generally, the curl of a vector) in a quantitatively correct manner requires the map factors.
The current implementation of
vorticityis correct if the underlying data is truly in Cartesian coordinates (even if unequally spaced). I suppose output from a cloud model would fit this description.However, when the data is not truly Cartesian (e.g., NWP output on some map projection), the third and fourth terms above are necessary for a quantitatively correct calculation. In many (most?) instances, the third and fourth terms will be quite small, certainly not noticeable qualitatively (say, when looking at contours), but one case where the difference is noticeable is when vorticity is computed from data on a lat/lon grid and examined in polar regions (where, due to convergence of the meridians, the dmx/dy term gets large).
For instance, compare the GEMPAK and MetPy-derived plots below (Python code is at https://gist.github.com/sgdecker/496f7ea7edd98b428ac3adab0b5e0842; GEMPAK is at https://gist.github.com/sgdecker/1a0bbf2e18dd0a2f9b4b02c83e07509b):
Look closely and you'll see, for instance, that the 4 contour to the left (south, technically) of the vort max in the Arctic Ocean differs noticeably.
A mathematical way to see this is to note that vorticity in spherical coordinates is given by (with theta as colatitude):
which expands to:
The final term (which corresponds to the dmx/dy term in GEMPAK's general expression) blows up at the pole, but if we are 1 degree away, the cotangent is order 10^2, so using U ~ 10 m/s, and r ~ 10^7 m, the term is on the order of 10^-4 s^(-1).
The divergence and laplacian operators will have the same issue.
dopplershift commentedon Jul 31, 2018
Thanks for the detailed analysis. It's what I figured, and definitely on our roadmap. Due to the nicely detailed analysis, let's leave this open for addressing the map factors, and we can close out #174.
(Putting the link to the NOAA writeup about map projections here: https://www.arl.noaa.gov/hysplit/cmapf-mapping-routines/)
[-]Generalize vorticity (and related functions) to allow for unequal grid spacing[/-][+]Generalize vorticity (and related functions) to allow for map projections[/+]deeplycloudy commentedon Jul 31, 2018
Nice resource in the ARL guide. I do want to suggest caution in retaining 3D coordinates as much as possible, as I advocated in the Cartopy PlateCaree saga - while the ARL guide has nice routines for 2D, and that's fine for maps, it can make forward/inverse and mixed-coordinate plotting difficult.
Anyway, proj.4 has a 3D-aware forward / inverse, and a quick search turns up partial derivative support through the proj_factors function - maybe that's useful in ensuring we account for what @sgdecker points out above. This issue on OSGeo also seems related.
sgdecker commentedon Aug 1, 2018
Since I mentioned divergence and the Laplacian operator, I went ahead and made plots of those as well to show the difference between GEMPAK and MetPy.
Here is 500-mb divergence:
Here is the Laplacian of 500-mb height:
Perusing the GEMPAK documentation, I found other calculations that also require taking derivatives of map factors to get absolutely correct results:
kgoebber commentedon Jun 12, 2019
Okay, so I've been looking into this and have the following to report on absolute vorticity...
Here is the difference between GEMPAK and MetPy v0.10; note the very small range of color values
The average difference is 0.0051 /s. If the following term is added
+ (uwnd / (6375000 * units.meter)) * np.tan(np.deg2rad(lats[:, None]))the difference is reduced by an order of magnitude to essentially zero.The average difference is -0.00042. Note: The single strip of color at the poles is because MetPy's derivative does not assess vorticity at the top and bottom of the array. Excluding those rows does not appreciably affect the mean difference calculations.
The additive factor depends on the u-component of the wind and latitude and it is not a universal or 3D solution to the issue as it will has the potential to be slightly different for some of the calculations (e.g., laplacian). Since the list is "relatively" small, this can be done for the various cases and added, however with our current implementation it would require the addition of potentially two arguments. One argument signaling to use (or not use) the spherical calculation and another that contains the latitude array. This would not be a problem for
absolute_vorticitysince it already requires the latitudes for the calculation ofcoriolis_parameter, therefore we would only need to add the spherical designation. I would additionally promote the default use of spherical to beTruesince that would be the much large use case for our current user base.sgdecker commentedon Jun 27, 2019
I haven't tried implementing this yet, but I think the general solution would be to accept the data CRS object as an optional argument, rather than more specific things such as latitude. If that argument is not provided, assume the grid is Cartesian, and set the map factors to one. Otherwise, use the CRS object to obtain the map factors via proj4 (as discussed by @deeplycloudy). In either case, use the formula in my original comment as the basis for the computation. The part that is unclear to me without actually trying this yet is how difficult obtaining the map factors from the CRS is.
jthielen commentedon Jun 28, 2019
This question is more related to the information needed for projection-aware deriviative-based calculations, rather than the implementation, but I still wanted to ask, since it may end up being useful for the implementation.
I've recently been planning out a collection of "dataset helpers" to flesh out any missing CF metadata in xarray Datasets, so that all of the MetPy-xarray integration functionality (projections, coordinates, units, and hopefully soon variable identification and the field solver) works as smoothly as possible. WRF output has been the motivating example (see #1004).
With that in mind, what would be the best set of information to standardize having in a DataArray and/or Dataset to ensure that these projection-aware calculations can be carried out? Is it sufficient to have
crs?or am I missing something else?
The goal, I would think, would be to be able to calculate something like vorticity with a call such as
and have everything else be obtained automatically from the DataArrays. This, if combined with some way of annotating the function and its arguments, would also really set the stage for the field solver.
dopplershift commentedon Jun 28, 2019
Correct, my ideal is to have the CRS information, as well as coordinates like x or lon, attached to the data arrays. If we can't get map factors from the CRS...that would be bad.
jthielen commentedon Jan 10, 2020
While I had hoped that #1260 (which allows for easy calculation of both nominal and actual grid deltas from a parsed DataArray) would be enough to make this feasible for quick implementation for v1.0, #1275 raised the point about grid-relative vs. earth-relative winds that I had not fully considered with respect to this issue. So,
With all this extra complexity, though, I unfortunately don't see a way that this can feasibly be resolved with a stable API for v1.0 in the immediate future. So, will this issue need to be punted from v1.0, and instead, for now, just use a simple xarray input/output wrapper that fills in the default grid spacing (
dxanddy) and coordinate-dependent parameters (fandlatitude), and calculates these kinematic functions as they are now?sgdecker commentedon Jan 10, 2020
The equation way back in the first comment is valid for grid-relative
uandv. If you used earth-relative winds on, say, a polar stereographic projection where North America is right-side up but Asia is upside down (so that grid-relative and earth-relative winds are nearly equal and opposite), I'd suspect adding the map factor terms would actually degrade the computation.38 remaining items
jthielen commentedon Feb 23, 2022
Not sure if this was the impression of other folks or not, but I previously conceptualized these map factor corrections as only occurring in derivatives on vector quantities. #2356 (comment) made me realize this does not seem to be the case; map factors need to be considered in any vector operation (so, also gradients of scalar fields, dot products, cross products, etc.). However, these corrections are much simpler (just coefficients on terms and they sometimes cancel out).
jthielen commentedon Feb 28, 2022
I keep on getting nerd sniped with this topic of "vector calculus on orthogonal coordinates," and while I don't have much new to show for it yet, one particular concern came up:
Do we need to delineate between vector component arguments that are scaled according to the covariant basis versus the normalized basis? (See https://en.wikipedia.org/wiki/Orthogonal_coordinates#Covariant_basis for discussion.) Is it safe to assume we will always have things in the normalized basis?
sgdecker commentedon Feb 28, 2022
For what it's worth, the Pielke textbook says (p. 129 of the second edition) "By convention, the equations [the set of equations solved in an NWP model] are written in the contravariant form, using the covariant differentiation operation..." but that is in the context of a whole bunch of derivations. I think in practice all of the vector quantities MetPy would be dealing with are normalized (e.g., if the u-component of the wind is given as 20 m/s, it really is 20 m/s at that point).
jthielen commentedon Mar 4, 2022
On the topic of frequently having this problem in the back of my mind, I ended up writing up a rough proposal for addressing this "grid-correct vector calculus" problem on the data model level: pydata/xarray#6331. I doubt that all would be practical on the timeline we want to have this functionality implemented in MetPy, but, it could be a cool future goal?
dopplershift commentedon Mar 4, 2022
👍 to all of that
dopplershift commentedon Dec 23, 2022
Implemented in #2743.