search

Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.21k stars 408 forks source link

Isallobaric Wind #641

Open jrleeman opened 6 years ago

jrleeman commented 6 years ago

AMS glossary cites Haurwitz, B. 1941. Dynamic Meteorology. 155–159., but I'm sure Holton or similar will have this as well.

ams2001glos-ie17

GEMPAK Docs

ISAL Isallobaric wind ISAL ( S ) = [ - DDT ( v (GEO(S)) ) / CORL, DDT ( u (GEO(S)) ) / CORL ]

jthielen commented 6 years ago

Does Bluestein 1992 or 1993 have this (since I think we have been citing those texts for many of the dynamics calculations), and if so, what page?

jrleeman commented 6 years ago

I wouldn't be surprised if it is in there and/or Holton. All of these are on the shelf in my office, feel free to borrow.

sgdecker commented 5 years ago

It's in Bluestein 1992 at the bottom of p. 173, the first term of (4.1.98).

sgdecker commented 2 years ago

This one is tricky because of the time derivative. By the way, it seems the current declarative interface is limited to one time (please correct me if I'm wrong). In GEMPAK, you would specify two times to plot this (the endpoints of the time interval). Similarly, in GEMPAK you would specify two levels for all the layer-based calculations (such as the thermal wind). Are there plans to generalize the declarative interface to allow these kinds of calculations, or should I open an issue?

dopplershift commented 2 years ago

@sgdecker I think we're seeing two options here in general:

  1. Do the calculation manually (maybe using one or more functions from metpy.calc) and hand the resulting DataArray to the declarative interface
  2. The solver (which I'm working on right now) will eventually gain some notion for levels/times. (Not sure how exactly that extends to the isallobaric case).

We're not necessarily wanting to have an entire domain specific language for specifying calculations like GEMPAK did. That's what Python is for.

sgdecker commented 2 years ago

Would possible implementations be limited to functions? I'm wondering about a calling convention like (my time strings may not be correct syntax): data.metpy.sel(time=slice('2021-08-24:12:00','2021-08-24:18:00')).ddt(dim='time') to compute the time derivative for every DataArray in the xarray Dataset data, or if data is a DataArray, just that DataArray.

I suppose this could be a function: time_der = ddt(data, slice('2021-08-24:12:00','2021-08-24:18:00'))

And the isallobaric wind (or any other computation that requires a time interval) would have a similar interface.

dopplershift commented 2 years ago

Couldn't ddt be done today (with the properly aggregated dataset) using metpy.calc.first_derivative(data_array, axis='time')?

sgdecker commented 2 years ago

Ahh, I looked at the GEMPAK Conversion Guide too quickly. That looks correct. BTW, what's the proper process to turn the "Tested against GEMPAK" column from red to green? Maybe open an issue that demonstrates the equivalence and a linked pull request that edits the table?

dopplershift commented 2 years ago

I'd say just the PR. If one is feeling particularly motivated, adding a test case for the relevant function that tests with values validated from GEMPAK? (Depends if that's meaningfully different from existing tests.)

sgdecker commented 2 years ago

Actually, metpy.calc.first_derivative isn't a drop in replacement for ddt, as first_derivative requires three points, but ddt uses only two: https://github.com/Unidata/gempak/blob/92f9a3a8ee667ec49a9082f44380e27f61ca716b/gempak/source/diaglib/df/dfddt.c

Given the coarse time resolution in most gridded datasets, this is somewhat of a limitation. For ddx and ddy, GEMPAK uses centered differences, so those cases are much closer to what first_derivative does. But ddt is a bit different.

sgdecker commented 2 years ago

It's not that one is right and the other is wrong, but they are different.

By the way, as I try to see how different the results are, I am noticing that first_derivative strips all the attributes from my DataArrary. The result then doesn't play nice with the declarative plotting interface. Here is what I am trying:

import datetime

import xarray as xr
from metpy.io import GempakGrid
from metpy.calc import first_derivative
from metpy.plots import ContourPlot, MapPanel, PanelContainer

MODEL = '/ldmdata/gempak/model/'
gem_data = GempakGrid(MODEL + 'nam/21082312_nam211.gem')

heights = gem_data.gdxarray(parameter='hght', level=500, coordinate='pres')
z_ds = xr.concat(heights, dim='time')

data = z_ds.sel(time=slice('2021-08-23-12:00','2021-08-24-12:00',2))
print(data)

dzdt = first_derivative(data, axis='time')
print(dzdt)

# Simple plot as sanity check works
#dzdt.sel(time='2021-08-24-00:00').plot()

# Try a more complete plot with declarative interface
ht_fall = dzdt.sel(time='2021-08-24-00:00').rename('dzdt') * 1e4

cntr = ContourPlot()
cntr.data = ht_fall
cntr.field = 'dzdt'
cntr.level = 500
cntr.time = datetime.datetime(2021,8,24)
cntr.contours = list(range(-20, 20, 4))
cntr.linecolor = 'black'
cntr.linestyle = 'solid'
cntr.clabels = True

panel = MapPanel()
panel.area = [-125, -74, 20, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = 'dzdt at {plot_time}'
panel.plots = [cntr]

pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]

pc.show()

with output:

<xarray.DataArray 'z' (time: 3, pres: 1, y: 65, x: 93)>
array([[[[5845.124 , 5847.4604, 5841.7964, ..., 5853.14  , 5853.14  ,
          5850.4204],
         [5846.116 , 5847.5884, 5847.668 , ..., 5856.372 , 5855.492 ,
          5855.4604],
         [5849.7964, 5850.036 , 5850.116 , ..., 5860.1323, 5859.284 ,
          5858.868 ],
         ...,
         [5617.7   , 5599.844 , 5584.1   , ..., 5662.724 , 5645.7163,
          5626.9644],
         [5608.0044, 5591.9243, 5578.532 , ..., 5653.3003, 5637.06  ,
          5619.188 ],
         [5601.3804, 5588.356 , 5578.804 , ..., 5645.076 , 5629.9243,
          5613.1084]]],

       [[[5839.3477, 5843.124 , 5843.204 , ..., 5851.876 , 5851.652 ,
          5855.876 ],
         [5844.308 , 5845.668 , 5845.796 , ..., 5859.796 , 5858.308 ,
          5856.9478],
         [5847.892 , 5848.6597, 5848.644 , ..., 5860.7397, 5861.3   ,
...
          5642.468 ],
         [5620.98  , 5615.204 , 5610.148 , ..., 5640.9478, 5637.028 ,
          5628.58  ],
         [5606.9   , 5604.4517, 5597.812 , ..., 5629.124 , 5624.9316,
          5616.772 ]]],

       [[[5843.7007, 5840.405 , 5842.6606, ..., 5852.7246, 5850.5327,
          5849.333 ],
         [5844.613 , 5843.509 , 5842.7246, ..., 5854.5166, 5854.837 ,
          5851.9565],
         [5845.2847, 5844.773 , 5843.253 , ..., 5857.189 , 5855.9565,
          5854.8047],
         ...,
         [5626.773 , 5617.3647, 5608.9326, ..., 5664.5166, 5658.341 ,
          5650.885 ],
         [5607.525 , 5598.693 , 5590.677 , ..., 5648.085 , 5641.6367,
          5633.381 ],
         [5587.8926, 5580.277 , 5573.5566, ..., 5632.1167, 5625.0767,
          5616.1006]]]], dtype=float32)
Coordinates:
  * time     (time) datetime64[ns] 2021-08-23T12:00:00 ... 2021-08-24T12:00:00
  * pres     (pres) int64 500
  * x        (x) float32 -4.226e+06 -4.145e+06 -4.064e+06 ... 3.17e+06 3.251e+06
  * y        (y) float32 2.035e+06 2.117e+06 2.198e+06 ... 7.155e+06 7.237e+06
Attributes: (12/17)
    crs_wkt:                        PROJCRS["unknown",BASEGEOGCRS["unknown",D...
    semi_major_axis:                6371200.0
    semi_minor_axis:                6371200.0
    inverse_flattening:             0.0
    reference_ellipsoid_name:       unknown
    longitude_of_prime_meridian:    0.0
    ...                             ...
    standard_parallel:              (25.0, 25.0)
    latitude_of_projection_origin:  0.0
    longitude_of_central_meridian:  -95.0
    false_easting:                  0.0
    false_northing:                 0.0
    grid_type:                      forecast
<xarray.DataArray (time: 3, pres: 1, y: 65, x: 93)>
<Quantity([[[[-2.50950566e-04 -1.19097674e-04  5.51689996e-05 ... -5.37165889e-05
    -3.87234158e-05  2.65158194e-04]
   [-6.63079156e-05 -4.16904026e-05 -2.94551143e-05 ...  1.79985894e-04
     1.37950756e-04  1.09411169e-04]
   [-3.59429253e-05 -2.80874747e-06  1.12802011e-05 ...  6.21880425e-05
     1.31830286e-04  1.44054272e-04]
   ...
   [ 7.13495325e-04  1.12608733e-03  1.51332714e-03 ... -2.71126076e-04
     7.97808612e-05  4.40894233e-04]
   [ 6.06271249e-04  9.99428078e-04  1.32312916e-03 ... -5.11514169e-04
    -5.44625741e-05  2.70549633e-04]
   [ 4.11642569e-04  8.38679561e-04  9.40726951e-04 ... -5.88531494e-04
    -1.75035265e-04  1.34978118e-04]]]

 [[[-1.64738408e-05 -8.16627785e-05  1.00029839e-05 ... -4.80934426e-06
    -3.01784939e-05 -1.25856753e-05]
   [-1.74006709e-05 -4.72174750e-05 -5.72148076e-05 ... -2.14753328e-05
    -7.58418330e-06 -4.05544705e-05]
   [-5.22189670e-05 -6.09164768e-05 -7.94361256e-05 ... -3.40666594e-05
...
     1.46117034e-04  2.76856599e-04]
   [-5.54967810e-06  7.83397533e-05  1.40561704e-04 ... -6.03626393e-05
     5.29706037e-05  1.64269341e-04]
   [-1.56108715e-04 -9.35081199e-05 -6.07356319e-05 ... -1.49993896e-04
    -5.61071325e-05  3.46317998e-05]]]

 [[[ 2.18002884e-04 -4.42278827e-05 -3.51630317e-05 ...  4.40979004e-05
    -2.16335720e-05 -2.90329545e-04]
   [ 3.15065737e-05 -5.27445475e-05 -8.49745009e-05 ... -2.22936560e-04
    -1.53119123e-04 -1.90520110e-04]
   [-6.84950087e-05 -1.19024206e-04 -1.70152452e-04 ... -1.30321361e-04
    -2.08858914e-04 -2.38116229e-04]
   ...
   [-5.03477874e-04 -7.20520020e-04 -9.38500298e-04 ...  3.12618679e-04
     2.12453206e-04  1.12818965e-04]
   [-6.17370605e-04 -8.42748571e-04 -1.04200575e-03 ...  3.90788891e-04
     1.60403781e-04  5.79890498e-05]
   [-7.23859999e-04 -1.02569580e-03 -1.06219822e-03 ...  2.88543701e-04
     6.28209997e-05 -6.57145182e-05]]]], '1 / second')>
Coordinates:
  * time     (time) datetime64[ns] 2021-08-23T12:00:00 ... 2021-08-24T12:00:00
  * pres     (pres) int64 500
  * x        (x) float32 -4.226e+06 -4.145e+06 -4.064e+06 ... 3.17e+06 3.251e+06
  * y        (y) float32 2.035e+06 2.117e+06 2.198e+06 ... 7.155e+06 7.237e+06
Traceback (most recent call last):
  File "/home/decker/classes/met433/new_labs/ddt_test.py", line 47, in <module>
    pc.show()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 589, in show
    self.draw()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 576, in draw
    panel.draw()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 846, in draw
    p.draw()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 1137, in draw
    self._build()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 1300, in _build
    x_like, y_like, imdata = self.plotdata
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 1121, in plotdata
    plot_x_dim = self.griddata.metpy.find_axis_number('x')
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 1086, in griddata
    data = self.data.metpy.parse_cf(self.field)
AttributeError: 'MetPyDataArrayAccessor' object has no attribute 'parse_cf'
dopplershift commented 2 years ago

Try not setting cntr.field since you're passing in a DataArray.

That's a good point about the 2 vs. 3 points. I'd be happy to merge something adding something functionally equivalent to ddt--naming it well will probably be the biggest challenge. The implementation (which can re-use @xarray_derivative_wrap and _process_deriv_args) is pretty much np.diff(f) / delta.

sgdecker commented 2 years ago

OK, I was overspecifying things. I also needed to drop the time specification, but my revised version still isn't happy:

import datetime

import xarray as xr
from metpy.io import GempakGrid
from metpy.calc import first_derivative
from metpy.plots import ContourPlot, MapPanel, PanelContainer
import cartopy.crs as ccrs

GRID211 = ccrs.LambertConformal(central_longitude=-95, standard_parallels=[25,25])
MODEL = '/ldmdata/gempak/model/'
gem_data = GempakGrid(MODEL + 'nam/21082312_nam211.gem')

heights = gem_data.gdxarray(parameter='hght', level=500, coordinate='pres')
z_ds = xr.concat(heights, dim='time')

data = z_ds.sel(time=slice('2021-08-23-12:00','2021-08-24-12:00',2))
dzdt = first_derivative(data, axis='time')

ht_fall = dzdt.sel(time='2021-08-24-00:00') * 1e4
ht_fall.attrs['crs'] = GRID211  # Error occurs with or without this line

cntr = ContourPlot()
cntr.data = ht_fall
cntr.level = 500
cntr.contours = list(range(-20, 20, 4))
cntr.linecolor = 'black'
cntr.linestyle = 'solid'
cntr.clabels = True

panel = MapPanel()
panel.area = [-125, -74, 20, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = 'dzdt at {plot_time}'
panel.plots = [cntr]

pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]

pc.show()

With output:

Traceback (most recent call last):
  File "/home/decker/classes/met433/new_labs/ddt_test.py", line 41, in <module>
    pc.show()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 589, in show
    self.draw()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 576, in draw
    panel.draw()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 846, in draw
    p.draw()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 1137, in draw
    self._build()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 1302, in _build
    kwargs = self.parent.plot_kwargs
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/plots/declarative.py", line 814, in plot_kwargs
    dataproj = self.plots[0].griddata.metpy.cartopy_crs
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/xarray.py", line 238, in cartopy_crs
    return self.crs.to_cartopy()
  File "/home/decker/local/miniconda3/envs/met433/lib/python3.9/site-packages/metpy/xarray.py", line 233, in crs
    raise AttributeError('crs attribute is not available.')
AttributeError: crs attribute is not available.

This is drifting away from the original issue, but I am attempting to understand in general how to combine calculations with GempakGrid objects and the declarative interface (if it is possible).

dopplershift commented 2 years ago

Change:

ht_fall.attrs['crs'] = GRID211

to

ht_fall.metpy.assign_crs(grid_mapping_name='lambert_conformal_conic',
                                           standard_parallel=[25, 25], longitude_of_central_meridian=-95)

See the xarray tutorial. This uses the names from CF metadata.

Should the parameters of the grid be within the GEMPAK file somewhere?

sgdecker commented 2 years ago

Thanks for the pointer to the documentation!

I think there are two issues: 1) The DataArrays from the GEMPAK reader do contain the key attributes (grid_mapping_name, standard_parallel, and longitude_of_central_meridian), but not in the form of a metpy_crs coordinate. I will draw upon the tutorial and your comment to make sure I can do a simple no-calculation example (just plot 500-mb heights) in the most straightforward way I can come up with. 2) first_derivative strips all those attributes away, so there is a need to re-attach them to the DataArray before trying to plot the result with the declarative interface.

sgdecker commented 2 years ago

I have success with the simple case!

import xarray as xr
from metpy.io import GempakGrid
from metpy.plots import ContourPlot, MapPanel, PanelContainer

MODEL = '/ldmdata/gempak/model/'
gem_data = GempakGrid(MODEL + 'nam/21082312_nam211.gem')

heights = gem_data.gdxarray(parameter='hght', level=500, coordinate='pres')
z_ds = xr.concat(heights, dim='time')

# This try fails with AttributeError: 'MetPyDataArrayAccessor' object has no attribute 'parse_cf'
# z_ds = z_ds.metpy.parse_cf()

z_ds = z_ds.metpy.assign_crs(grid_mapping_name=z_ds.grid_mapping_name,
                             standard_parallel=z_ds.standard_parallel,
                             longitude_of_central_meridian=z_ds.longitude_of_central_meridian,
                             latitude_of_projection_origin=z_ds.latitude_of_projection_origin)

z500 = z_ds.metpy.sel(time='2021-08-23 12:00', vertical=500)

cntr = ContourPlot()
cntr.data = z500
cntr.contours = list(range(5400, 6000, 60))
cntr.linecolor = 'black'
cntr.linestyle = 'solid'
cntr.clabels = True

panel = MapPanel()
panel.area = [-125, -74, 20, 55]
panel.projection = 'lcc'
panel.layers = ['states', 'coastline', 'borders']
panel.title = '500-mb heights'
panel.plots = [cntr]

pc = PanelContainer()
pc.size = (15, 15)
pc.panels = [panel]

pc.show()

But I am confused why parse_cf doesn't work. Using assign_crs seems like something gdxarray should do. Am I overlooking anything here?

I will move on to the ddt/first_derivative case next.

dopplershift commented 2 years ago

@sgdecker first_derivative, like pretty much all operations, calculations, math, etc. done on DataArray instances dumps attributes (no easy way to prevent this or to propagate them correctly). This is why parse_cf converts that metadata into an object that gets stuffed into the coordinates.

The reason parse_cf doesn't work is because that only exists on Datasets, because of the way the grid metadata is specified in the CF specification. For CF, the grid mapping parameters (e.g. standard_parallel, earth_radius) are attributes of a scalar variable, e.g. MyProjection, then all of the data variables reference that by having an attribute grid_mapping_name with a value of "MyProjection". I agree that assign_crs should be called by all of the xarray functions for the GEMPAK code--would you open a separate issue about that?

jthielen commented 2 years ago

Hope you both don't mind me chiming in here:

That's a good point about the 2 vs. 3 points. I'd be happy to merge something adding something functionally equivalent to ddt--naming it well will probably be the biggest challenge. The implementation (which can re-use @xarray_derivative_wrap and _process_deriv_args) is pretty much np.diff(f) / delta.

This won't be able to use @xarray_derivative_wrap as-is, since the output won't be on the same points/coordinate as the input (size will be one smaller, so have to decide if this acts as forward, backward, or centered difference and determine the new coordinate array accordingly).

Also, not to prematurely generalize, but this strikes me as something that falls under the alternative calculations (https://github.com/Unidata/MetPy/issues/142) area. Perhaps some careful thought here in the API can lead the way for other alternative calculations?

  1. The DataArrays from the GEMPAK reader do contain the key attributes (grid_mapping_name, standard_parallel, and longitude_of_central_meridian), but not in the form of a metpy_crs coordinate. I will draw upon the tutorial and your comment to make sure I can do a simple no-calculation example (just plot 500-mb heights) in the most straightforward way I can come up with.

I'd consider this a bug in the GEMPAK reader. I'm not aware of attributes containing the grid mapping details being valid on the variable (DataArray) itself. Instead, since we're building this from a pyproj.CRS anyway, we should just include it in the form of a metpy_crs. I'll open another issue for this.

  1. first_derivative strips all those attributes away, so there is a need to re-attach them to the DataArray before trying to plot the result with the declarative interface.

Correct, calculations stripping attributes is intended behavior, because it is not safe to assume attributes are still valid on the output (they represent different physical variables). At least in my understanding, things that are preserved across operations are necessarily more coordinate-like, so it makes sense (like metpy_crs) to make them actual coordinates.

But I am confused why parse_cf doesn't work. Using assign_crs seems like something gdxarray should do. Am I overlooking anything here?

Since CF grid mappings work by reference to another variable on the dataset, it doesn't make sense for parse_cf to exist for DataArrays (since the reference to the defining attributes is lost). As mentioned above, creating the metpy_crs is indeed something the GEMPAK reader should do internally.

jthielen commented 2 years ago

Oh wait, it looked like @dopplershift beat me to it on those points.

sgdecker commented 2 years ago

OK, the first_derivative calculation is fine, too, since the metpy_crs coordinate is retained even though the attributes go away (which is the point of metpy_crs at least in part).

Actually, the comparison is a lot closer than I thought it would be (at least for this case): ddt ddt

The one use case where it would still be worthwhile to have a ddt equivalent would be if you didn't have any time points to spare (for instance, you have model output at 3-hr intervals and want to compute ddt (or something based on ddt like the isallobaric wind) over a 3-hr interval).

And I see the new comments as well. I think the action items are:

  1. Update the GEMPAK Comparison Guide to state that first_derivative is almost a ddt replacement
  2. Implement ddt
  3. Implement an isallobaric wind function
  4. Open an issue to add the creation of metpy_crs to ... well, Jon just took of this one!
spacekace commented 1 month ago

Hi all,

I wanted to follow up on the ddt and isallobaric wind function action items from August 2021. I am looking to plot the various components of the ageostrophic wind including the isallobaric/isallohypsic wind. There is already an existing function for the inertial advective wind which is great - thank you!

I have some preliminary working code for both the isallobaric and inertial diabatic terms, but upon looking through this discussion on ddt I am unsure of how accurate some of my calculations are. If ddt and the isallobaric wind have already been developed, I would really appreciate using those functions as opposed to reinventing the wheel.

At this point, I wonder if it would be easier to simply calculate both the isallobaric and inertial advective terms and subtract those from the ageostrophic wind to get the inertial diabatic component since calculating the diabatic heating rate is quite a process.