How to use calculate shearing & stretching deformation?

[kminzoo]

Hi! I'm trying to calculate shearing & stretching deformation by using the MetPy functions, mpcalc.stretching_deformation() and mpcalc.shearing_deformation().

but i'm not sure that i correctly use them, because the result was different from the result of another method to calculate deformations...

the code is here.

import metpy
import metpy.calc as mpcalc

u = xr.open_dataset("u850_korea_202008.nc").squeeze()
v = xr.open_dataset("v850_korea_202008.nc").squeeze()

u_name= 'u'
v_name= 'v'
vo_name='vo'
lon_name='longitude'
lat_name='latitude'

wind = xr.merge([u, v])
wind = wind.rename({u_name: "u", v_name: "v", lat_name: "lat", lon_name: "lon"})

E = mpcalc.stretching_deformation(wind.u, wind.v)
F = mpcalc.shearing_deformation(wind.u, wind.v)

the dataset 'wind' looks like this.

[dopplershift]

That code using the functions looks correct. What method were you using for comparison? What version of MetPy are you using?

[kminzoo]

The version of MetPy is 1.4.
I use functions of a 'DynamicoPy' package as the another method.

The code is here.

import xarray as xr

u = xr.open_dataset("u850_korea_202008.nc").squeeze()
v = xr.open_dataset("v850_korea_202008.nc").squeeze()

u_name= 'u'
v_name= 'v'
vo_name='vo'
lon_name='longitude'
lat_name='latitude'

wind = xr.merge([u, v])
wind = wind.rename({u_name: "u", v_name: "v", lon_name: "lon", lat_name: "lat"})

#This part is different from the first method.
E = compute_stretching_xr(wind)
F = compute_shearing_xr(wind)

And the functions is here.

def get_dx_dy(lon, lat):

   dlon = lon[1] - lon[0]  # resolution in longitude in deg
   dlat = lat[1] - lat[0]  # resolution in latitude in deg
   R = 6371000  # Earth radius
   dy = R * dlat * np.pi / 180  # Vertical size of a cell
   # Horizontal size of a cell depending on latitude
   lat_rad = lat * np.pi / 180  # Latitudes in rad
   r = np.sin(np.pi / 2 - abs(lat_rad)) * R
   dx = r * dlon * np.pi / 180
   dx = np.transpose([(dx[1:] + dx[:-1]) / 2] * (len(lon) - 1))
   return dx, dy

def compute_stretching_xr(xarray, u_name="u", v_name="v", lon_name="lon", lat_name="lat"):

   xarray = xarray.rename(
       {lat_name: "lat", lon_name: "lon", u_name: "u", v_name: "v"}
   ).squeeze()

   lon, lat = xarray.lon.values, xarray.lat.values
   dx, dy = get_dx_dy(xarray.lon.values, xarray.lat.values)

   u = xarray.u
   v = xarray.v
   E = (
       u.sel(lat=lat[:-1], lon=lon[1:]).values
       - u.sel(lat=lat[:-1], lon=lon[:-1]).values
   ) * 1 / dx[:] - (
       v.sel(lat=lat[1:], lon=lon[:-1]).values - v.sel(lat=lat[:-1], lon=lon[:-1])
   ) * 1 / dy
   E = xr.DataArray(E).interp_like(xarray, kwargs={"fill_value": "extrapolate"})
   E.attrs["units"] = "s-1"
   E = E.rename({"lat": lat_name, "lon": lon_name})
   return E

def compute_shearing_xr(xarray, u_name="u", v_name="v", lon_name="lon", lat_name="lat"):

   xarray = xarray.rename(
       {lat_name: "lat", lon_name: "lon", u_name: "u", v_name: "v"}
   ).squeeze()

   lon, lat = xarray.lon.values, xarray.lat.values
   dx, dy = get_dx_dy(xarray.lon.values, xarray.lat.values)

   u = xarray.u
   v = xarray.v
   F = (
       v.sel(lat=lat[:-1], lon=lon[1:]).values
       - u.sel(lat=lat[:-1], lon=lon[:-1]).values
   ) * 1 / dx[:] - (
       u.sel(lat=lat[1:], lon=lon[:-1]).values - u.sel(lat=lat[:-1], lon=lon[:-1])
   ) * 1 / dy
   F = xr.DataArray(F).interp_like(xarray, kwargs={"fill_value": "extrapolate"})
   F.attrs["units"] = "s-1"
   F = F.rename({"lat": lat_name, "lon": lon_name})
   return F

[dopplershift]

So those implementations look like they use:

1. Grid spacing calculated assuming a spherical Earth
2. Uses a 2-point finite differencing scheme

(1) is similar to how MetPy worked pre-1.4, though MetPy by default assumes a WGS-84 spheroid. Regarding (2), MetPy uses a 3-point (centered) scheme that accounts for non-uniform grid spacing.

In MetPy 1.4, we properly calculate all the spatial derivative terms in orthogonal curvilinear coordinates (essentially accounting for the fact that the unit vectors change across the domain). See #1963. This is more correct when performing spatial derivatives on a spherical Earth.

How big are the differences you're observing? You might want to compare against MetPy <1.4 just to see how that compares with your expectation. Regardless, there are enough differences between the two approaches that I would only expect qualitative agreement, not numerically exact agreement.