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.24k stars 413 forks source link

Add capability for log _interp to interpolate to different values on multidimensional arrays. #524

Open mwilson14 opened 7 years ago

mwilson14 commented 7 years ago

Right now, log_interp can only interpolate to a multidimensional array of values if it's trying to interpolate all of them to the same value (say, sigma-to-pressure interpolation of model output to a 200mb surface). It would be nice to be able to interpolate multi-dimensional data to a 'surface' of different values (for example, interpolating a 2d gridded array of lcl pressures output by the lcl function to height using a 3d array of heights).

dopplershift commented 7 years ago

That sounds....incredibly tricky to get the numpy indexing right, though useful.

jthielen commented 6 years ago

With DataArray.interp(), I was able to do this with non-logarithmic interpolation (so long as we don't have a time coordinate, which is broken at the moment but should be fixed in v0.10.8). Here's a quick example (that could probably be cleaned up a little):

temp = data.isel(isobaric3=-1, time1=0)['temperature']
rh = data.isel(isobaric3=-1, time1=0)['relative_humidity']
height = data['height'].isel(time1=0)
dewpoint = mpcalc.dewpoint_rh(temp, rh / 100)

lcl_pres, lcl_temp = mpcalc.lcl(
    1000 * np.ones_like(dewpoint) * units.hPa,
    temp,
    dewpoint
)
lcl_pres = xr.DataArray(
    lcl_pres,
    coords=temp.coords,
    dims=temp.dims,
    name='lcl_pres',
    attrs={'units': lcl_pres.units}
)
lcl_height = height.interp(
    isobaric3=lcl_pres,
    latitude=lcl_pres['latitude'],
    longitude=lcl_pres['longitude']
)
print(lcl_height)
<xarray.DataArray 'height' (latitude: 81, longitude: 131)>
array([[2366.745918, 2157.203547, 1838.610165, ...,  220.58734 ,  249.562333,
         291.711661],
       [1955.285252, 1946.359996, 1781.497647, ...,  206.770545,  222.959318,
         243.481029],
       [1986.662478, 1893.573295, 1811.840817, ...,  190.616433,  202.580812,
         231.03481 ],
       ...,
       [ 539.128245,  529.080681,  519.104705, ...,  780.501756,  796.804493,
         796.198357],
       [ 520.75983 ,  516.120474,  503.938927, ...,  797.27822 ,  742.306889,
         761.436377],
       [ 508.236845,  498.527472,  511.411508, ...,  890.443381,  795.858029,
         749.889825]])
Coordinates:
    time1      datetime64[ns] 2017-09-05T12:00:00
    reftime    datetime64[ns] 2017-09-05T12:00:00
    crs        object Projection: latitude_longitude
    isobaric3  (latitude, longitude) float32 774.051 793.80334 825.368 ...
  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 48.0 47.5 47.0 46.5 ...
  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 252.0 252.5 253.0 ...
Attributes:
    long_name:                      Geopotential height @ Isobaric surface
    units:                          m
    Grib_Variable_Id:               VAR_0-3-5_L100
    Grib2_Parameter:                [0 3 5]
    Grib2_Parameter_Discipline:     Meteorological products
    Grib2_Parameter_Category:       Mass
    Grib2_Parameter_Name:           Geopotential height
    Grib2_Level_Type:               100
    Grib2_Level_Desc:               Isobaric surface
    Grib2_Generating_Process_Type:  Forecast

This should be easy to wrap for log interpolation as well.