Replies: 4 comments 3 replies
-
|
This sounds really exciting! Many Xarray users, including myself, would find this useful. I'd love to see a prototype -- have you thought about what the implementation could look like yet? |
Beta Was this translation helpful? Give feedback.
-
|
I had some very preliminary thoughts on implementation, such as
but haven't gotten around to any sort of prototype due to
So, hopefully with some discussion/feedback here I can work up to such a prototype, since #5692 looks really close from my outside point-of-view! |
Beta Was this translation helpful? Give feedback.
-
|
I'll be honest, as a programmer who hasn't done anything with vectors beyond converting wind u/v components to wind direction and speed, a lot of this proposal is over my head. That said, if we're talking about CRS, pyproj, and rioxarray then @snowman2 is another person who has be brought into the conversation. Regarding geoxarray, I wanted geoxarray to serve two purposes:
A lot of the above is based on my limited 2D single-coordinate-per-pixel satellite data experience. If geoxarray needs to be limited to "only" purpose 1 and maybe the pyproj transformation in 2 and resampling gets moved to another package or absorbed into rioxarray that is fine by me, especially if it serves the community better. I say this because it wasn't clear to me but it sounds like what you're proposing would require at least 2 new libraries: the geolocation/geoxarray-style and the vector stuff with some of those being built on or being part of xarray core. Does that sound right? |
Beta Was this translation helpful? Give feedback.
-
|
Thank you for the elaboration on the geoxarray side of things!
I was thinking just a single package, but you might be on to something about needing multiple libraries to accomplish these goals. The reason why I envisioned just a single package is the central role that
At this very early stage, I don't think a firm decision needs to be made either way! I could prototype either approach once I get around to prototyping. But, if this "xfield" idea gains traction, then we should certainly coordinate among interested packages to have an agreed upon place in the data model where the CRS is stored and interfaced with, preferably in one central (i.e., lightweight) package's implementation of a custom index. |
Beta Was this translation helpful? Give feedback.
-
|
@jthielen This sounds indeed very interesting. Our current approach is indeed somewhat hacky at the moment. From my xgcm perspective I am unsure about
Does this mean the DataVectorField would 'know' how to compute certain vector calculus operators? I haven't completely absorbed this logic here in detail, but would love to talk more. I am also curious as to what @TomNicholas and @rabernat think? |
Beta Was this translation helpful? Give feedback.
-
|
As questions have been asked, I've had some subtle changes of thought, so I apologize for this rambling mini-essay of a response!
Yes, so long as there are true orthogonal coordinates, that was the intent! Here's some pseudo-code/incomplete Python for what I'd imagine on that end: class DataVectorField:
def divergence(self):
return sum(
partial_derivative(
(
self.scale_factor_jacobian_determinant
* self.get_map_factor(field_dim)
* self[field_dim]
),
field_dim
)
for field_dim in self.field_dims
) / self.scale_factor_jacobian_determinant
@verify_dimensionality(2)
def k_curl(self):
return (
partial_derivative(
self.get_map_factor(self.field_dims[1]) * self[self.field_dims[1]],
self.field_dims[0]
)
- partial_derivative(
self.get_map_factor(self.field_dims[0]) * self[self.field_dims[0]],
self.field_dims[1]
)
) / self.scale_factor_jacobian_determinant
@verify_dimensionality(3)
def curl(self):
return sum(
(
self.get_map_factor(self.field_dims[i])
* self.get_map_factor(self.field_dims[j])
* (
partial_derivative(
self[self.field_dims[j]] / self.get_map_factor(self.field_dims[j]),
self.field_dims[i]
)
- partial_derivative(
self[self.field_dims[i]] / self.get_map_factor(self.field_dims[i]),
self.field_dims[j]
)
)
* self.normalized_basis[k]
)
for k, j, i in ((0, 2, 1), (1, 0, 2), (2, 1, 0))
)
def get_map_factors(self, field_dim):
if f'map_factors_{field_dim}' in self.coords:
# Cached
# This would ideally use metadata (i.e., cf-xarray), but here name gets
# the point across
return self.coords[f'map_factors_{field_dim}']
elif self.field_index.has_pyproj_crs() and self.is_horizontal(field_dim):
# Missing, but we can compute
self.cache_map_factors_from_pyproj()
return self.coords[f'map_factors_{field_dim}']
else:
# No scale factors available
raise ValueError('Cannot determine map factors.')
@property
def scale_factor_jacobian_determinant(self):
return np.stack(
[self.get_map_factors(field_dim)**-1 for field_dim in self.field_dims],
axis=0
).prod(axis=0)
@requires('pyproj')
def cache_map_factors_from_pyproj(self):
dims = self.field_index.horizontal_geographic_dims()
meridional_dim, parallel_dim = dims
factors = pyproj.Proj(self.field_index.pyproj_crs).get_factors(
self.get_auxilliary_coordinate('longitude', dims=dims),
self.get_auxilliary_coordinate('latitude', dims=dims)
)
self.coords[f'map_factors_{meridional_dim}'] = xr.Variable(
dims=dims,
data=factors.parallel_scale
)
self.coords[f'map_factors_{meridional_dim}'] = xr.Variable(
dims=dims,
data=factors.parallel_scale
)
# other things forward to the field_index, which is a bit "magic" at this
# point...I have little idea of what the implementation would look like,
# just what it does/provides (wraps other xarray indexes, tracks field
# dimensions, maintains a mapping between array dimensions and field dimensions,
# and if geographic, stores a pyproj.CRS and labels which field dimensions have
# "X" and "Y" semantics)
# `normalized_basis` is a tuple of special `BasisVector` objects which turn
# scalar fields into a vector whose single non-empty/non-zero component
# is in that basis direction.
def partial_derivative(field, axis):
# For non-periodic, unstaggered grids, can just vendor MetPy's implementation
# of Bowen et al. (2005), which handles variable grid spacing, and is a
# three-point centered difference on interior points and three-point one-sided
# difference on edges
# https://github.com/Unidata/MetPy/blob/v1.2.0/src/metpy/calc/tools.py#L950-L1041
return metpy.calc.first_derivative(field, axis=axis)
# Periodic, unstaggered grids can use a simplification of the previous with
# cyclic indexing
# This technically also works on staggered grids, but gives local quantities
# rather than cell-center/cell-corner from cell-wall like most finite-volume-type
# derivatives in models behave. So, we get to a really subtle point of "what
# kind of derivative do we want?"
Understandable--especially since I haven't formulated all the logic yet either! How something like this would play along with xgcm is a particular uncertainty for me, since
So, after thinking about it some more, how I envisioned this "xfield" idea sits in the middle of different levels of xgcm functionality, which makes it a bit of a sticky problem. Also, in my desire for a slick API, I realized I totally glossed over the subtleties of "what kind of derivative", which is a crucial point. With all that, perhaps it would be better to split the idea up a bit more? Just brainstorming in the open with this:
These levels could also guide implementation bottom-up (e.g., provide the data model for what xgcm needs now and what's envisioned for other levels first, then worry about the vector calculus primitives, then worry about the nice API, and so on)? |
Beta Was this translation helpful? Give feedback.
-
|
I am trying to implement some atmospheric sciences physics formulas and I just hit a weird quirk with my current workflow (using metpy to calculate divergence and vorticity of an Xarray.Dataset and then using SHTools to calculate the Spherical Harmonics). I would love to contribute a bunch of calculation functions to Metpy or Xarray but I am not sure how to handle basic partial derivatives on Xarray Vector fields (say ds.UU with coords lon, lat, pres and time) or more complicated ones with 2 dimensional coords (lon=(rlon,rlat)). Anyone has any updates on this? I would not mind writing a documentation page for such use cases 🙂 |
Beta Was this translation helpful? Give feedback.
-
Related to an earlier "ideas" discussion I posted on vectors in xarray, I've spent some time trying to brainstorm improved ways of representing scalar and vector fields (like those frequently encountered in models in the geosciences) and enabling easy vector calculus operations with them. I wanted to share some of those thoughts to see if such functionality already exists in the broader xarray ecosystem that I'm not yet aware of, and if not, if there is community interest for the kind of extension package I had in mind.
Tagging some folks who may have thoughts on this based on prior related discussions and overlapping tools/libraries: @keewis @rabernat @michaelavs @erogluorhan @benbovy @willirath @dcherian @eteq @djhoese @jbusecke @TomNicholas @dopplershift
Also, if there is a better forum for this type of brainstorming/proposal over xarray's Discussions, please let me know!
The Problem
Overview
xarray.DataArrayprovides a wonderful data model for labeling dimensions of an n-dimensional array and providing coordinates/indexes associated with those dimensions. However, when working with scalar and vector fields, dimensions no longer have just one semantic meaning--instead, we haveWhen working with scalars alone (e.g., like those considered in the CDM or CF Conventions), this distinction is often unnecessary at the data model level. However, it becomes critical when representing vector fields and/or performing vector calculus in non-Cartesian coordinates (like those that occur with nearly all geospatial data, such as rectilinear and projected grids). In the latter case, the coordinates associated with the field dimensions carry with them basis vectors and scale factors that would ideally have a high-level representation in code.
Motivating Example
Given an
xarray.Datasetcontaining model data over four dimensions (time, vertical, y, and x) with both scalars and vector field components, how can we represent the operations that make up something like vertical vorticity advection?With existing xarray functionality, a package could definitely take the approach of bundling all such calculations into their own functions in terms of scalars and vector components,
however, this is rather inflexible. Beyond not representing vectors as vectors (only components, which could be subtly error-prone), it relies upon the library developers to enumerate out all such fundamental calculations that users may wish to have. Instead, I'd love to see the operations themselves be able to be represented in code, so something like:
To do this, we would need:
DataArrayto make a "Scalar Field" type (I think a customFieldIndexand axfieldaccessor should do the trick!)Envisioned Use Cases
MetPy
A long requested enhancement to MetPy's suite of atmospheric kinematics calculations is the robust handling of projection scale factors. Presently (v1.2 and earlier), MetPy's function signatures look a bit like
and are designed to flexibly handle either
pint.Quantitys orxarray.DataArrays (in the latter case, optional arguments are automatically determined from coordinate information). However, to handle scale factors in the same way (while also making it as easy as possible for users to provide information they already have rather requiring prerequisite boilerplate code), this "blows up" to something likeBeyond that, the underlying conditional inclusion of scale factors throughout the calculation stack significantly increases the code complexity and verbosity of internal implementations. If the vector calculus operations and data representations themselves could be effectively encapsulated, this would be much easier to implement, test, maintain, and extend.
xgcm
xgcm does a lot of great things with scalar and vector handling on staggered grids, however, it seems to suffer due to the lack of a true vector data representation (e.g., xgcm/xgcm#410). If a
xfieldpackage containing aDataVectorFieldtype were to exist, I'd think much of xgcm's API could be cleaned-up/streamlined. What would such aDataVectorFieldlook like for a staggered grid? This is the kind of repr I'd imagine for a 2D vector on a C-grid implicitly in 3D space (like the common case of horizontal wind in the atmosphere):However, in memory, it would look a lot like a
Dataset(or aDataArraygroup-by along a stack dimension, as in #5775), but withFieldIndexPandasIndexand the likely futurexgcmIndexto support true indexing operations, rather than just field alignment and scaling propertiesdata_vars/componentsare aligned/expressed in terms of the axes of the field coordinate system (stored in theFieldIndex)geoxarray /
CRSIndex/WCSIndexThis project would also aim to handle the "coordinate system representation problem" expressed in #3620 and related packages such as geoxarray, and potentially be able to subsume some of their functionality depending on the approach. However, I don't think this
xfieldidea would fit itself within the existinggeoxarray, because this is meant to also work on scalar and vector fields in a more general sense. Neither wouldgeoxarraybe obviated, since it has intended functionality outside of just "carrying the CRS around," such as format/representation conversions.xoak
If xoak (or similar package) were to implement the long-coveted "2D longitude/latitude"
Index, it would be natural to wrap that with the coordinate system handling provided here. For example, a horizontal wind vector with west-east and south-north components, but with data existing on an unspecified projected grid, could look likeLimits of Scope (at least to begin)
Tensor Fields
If this handles scalar fields and vector fields, would it not make sense to generalize to tensor fields? With popularity of things like
einops, there's likely a place for it within the community.That being said, vectors map to lists, so it's easy to conceptualize as a data structure, and in that vein, keeping array dimension <-> field coordinate <-> vector component associations is straightforward. That no longer holds with tensors in general (especially mixed tensors), so that would be a whole other level of abstraction to handle. Also, tensors in curvilinear coordinates is way beyond my areas of expertise, so I can't even conceptualize what this would look like in a user-friendly way, which puts this functionality beyond practicality unless someone else with that expertise takes the lead on it.
Unstructured Grids (uxarray)
There is definitely some overlap with the "Application of calculus operations, including divergence, curl, Laplacian and gradient" stated goal. However, unless we can find a suficiently robust way to arbitrarily relate
then I don't think the approach used here could work. The "differential operations along orthogonal coordinates" is kind of central to the outlined approach to vector calculus here, and trying to generalize that to differential operators on meshes (where discritization is often based on cell geometry and integral-limit definitions of operators, rather than building up from differences in orthogonal directions) is a much more difficult problem. That being said, I could see collaboration to establish harmonious APIs to be a good goal! And perhaps there is something I'm missing that would allow unstructured grid handling here as well.
Integral Calculus, not just Differential
This really should be within scope for this kind of package, particularly with integration with
xgcm. However, my mesoscale meteorology bias shines through in that I don't really use much in the way of integration-type operations, so it's hard to motivate (or test from experience) such functionality from my own use cases.This is definitely an area (pun intended) where it'd be great to have community engagement to develop.
CRS Format Conversions, Interpolations, Re-Projection, and other GIS-type operations
This seems like a better fit for
geoxarray,xesmf,rioxarray, or other community packages.With CRS format/storage, this package should either pick a single, consensus-based CRS representation to use internally, or be generic. Given the need to "automagically" compute scale factors, it almost certainly is going to be the former, and that representation is going to be
pyproj.CRS.Other GIS operations would definitely bloat the scope of this package too much. I think here we'd just need to be careful to transparently pass through everything that the indexes wrapped within the
FieldIndexneed to work properly, and then let those other packages do all the work.Beta Was this translation helpful? Give feedback.
All reactions