Kinematic calculations on lat-lon grid

[mccrayc]

Currently, kinematic calculation functions (vorticity, advection, etc.) take arguments of dx, dy grid spacing in distance units. I'd like to calculate these functions for my data on a 0.5°x0.5° lat-lon grid. It would be useful to be able to do these calculations with dx, dy input as lat/lon degrees.

I.e. v_vorticity(u, v, 0.5,0.5, latlon=Y)

[dopplershift]

I think my preference would be to rely on our units infrastructure to look at the units on dx, dy, rather than a boolean flag.

Would also be nice to take the the full grids--though really all of the stuff involving derivatives and knowing which axis is which dimension should really be using xarray. Always more to be done. 😁

[mccrayc]

Ahh yeah - definitely agree on using units! Perhaps NCL docs for relative vorticity could be useful inspiration. I think allowing for calculations over lat/lon grids could go a long way to providing some more GEMPAK-esque functionality.

[dopplershift]

If you come up with some way to make it work, please let us know--we'd love to have a Pull Request. 😄

[kgoebber]

Here is some code to calculate the distance between latitude and longitude points, provided it is on a "clean" grid, a la GFS data. The example below makes some fake 0.5 degree lat/lons (similar to what you get from the 0.5 GFS) to calculate the distance in meters along the x and y (lon/lat) dimensions. It gives back two 2D arrays for easy calculations for a derivative over a grid space and appears to work well with the numpy gradient function.

import numpy as np

def calc_dx_dy(longitude,latitude):
    ''' This definition calculates the distance between grid points that are in
        a latitude/longitude format.
        
        Equations from:
        http://andrew.hedges.name/experiments/haversine/

        dy should be close to 55600 m
        dx at pole should be 0 m
        dx at equator should be close to 55600 m
        
        Accepts, 1D arrays for latitude and longitude
        
        Returns: dx, dy; 2D arrays of distances between grid points 
                                    in the x and y direction in meters 
    '''
    dlat = np.abs(lat[1]-lat[0])*np.pi/180
    dy = 2*(np.arctan2(np.sqrt((np.sin(dlat/2))**2),np.sqrt(1-(np.sin(dlat/2))**2)))*6371000
    dy = np.ones((latitude.shape[0],longitude.shape[0]))*dy

    dx = np.empty((latitude.shape))
    dlon = np.abs(longitude[1] - longitude[0])*np.pi/180
    for i in range(latitude.shape[0]):
        a = (np.cos(latitude[i]*np.pi/180)*np.cos(latitude[i]*np.pi/180)*np.sin(dlon/2))**2
        c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a) )
        dx[i] = c * 6371000
    dx = np.repeat(dx[:,np.newaxis],longitude.shape,axis=1)
    return dx, dy

lat = np.arange(90,-0.1,-0.5)
lon = np.arange(0,360.1,0.5)
print(lat)
print(lon)

dx, dy = calc_dx_dy(lon,lat)

print(dx)
print(dy)

[dopplershift]

At the risk of opening up a can of worms: what about other projections? Just project with cartopy?

[kgoebber]

Just haven't tried to tackle that yet. Let me see what I can come up with if I use some NAM like data. My guess is that as long as you can get lat/lon values, I'll be able to abstract the code to work for any 2D lat/lon grid.

[kgoebber]

Here is modified code to accept 2D lat/lon arrays, could easily add check for 1D arrays and make them 2D within definition. It gives the same output for the 1D arrays in the previous example code above. To test non-GFS like data I used the NAM lat/lon arrays pulled through Siphon. After emailing with Sean, he pointed me to the online documentation for the different grids. I have been able to confirm that the code is giving the correct distances at the different latitudes. The NAM grid spacing (for the 218 grid) should be 12 km at 35N, while at 25N it is 12.19 km.

http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID218

I'm sure the code could be sped up from how it is written, but I believe this would get the job done for calculating the distance between grid points for 2D lat/lon grids.

http://nbviewer.jupyter.org/url/fujita.valpo.edu/%7Ekgoebber/NAM_distance_lat_lon.ipynb

def calc_dx_dy(longitude,latitude):
    ''' This definition calculates the distance between grid points that are in
        a latitude/longitude format.
        
        Equations from:
        http://andrew.hedges.name/experiments/haversine/
        
        Accepts, 1D arrays for latitude and longitude
        
        Returns: dx, dy; 2D arrays of distances between grid points in the x and y direction in meters 
    '''
    dx = np.empty(latitude.shape)
    dy = np.zeros(longitude.shape)
    
    dlat_x = (latitude[1:,:]-latitude[:-1,:])*np.pi/180
    dlon_x = (longitude[1:,:]-longitude[:-1,:])*np.pi/180
    
    dlat_y = (latitude[:,1:]-latitude[:,:-1])*np.pi/180
    dlon_y = (longitude[:,1:]-longitude[:,:-1])*np.pi/180
    
    #print(dlat_y.shape)
    #print(dlon_y.shape)
    for i in range(latitude.shape[1]):
        for j in range(latitude.shape[0]-1):
            a_x = (np.sin(dlat_x[j,i]/2))**2 + np.cos(latitude[j,i]*np.pi/180) * np.cos(latitude[j+1,i]*np.pi/180) * (np.sin(dlon_x[j,i]/2))**2
            c_x = 2 * np.arctan2(np.sqrt(a_x), np.sqrt(1-a_x))
            dx[j,i] = c_x * 6371229
    dx[j+1,:] = dx[j,:]
    
    for i in range(latitude.shape[1]-1):
        for j in range(latitude.shape[0]):       
            a_y = (np.sin(dlat_y[j,i]/2))**2 + np.cos(latitude[j,i]*np.pi/180) * np.cos(latitude[j,i+1]*np.pi/180) * (np.sin(dlon_y[j,i]/2))**2
            c_y = 2 * np.arctan2(np.sqrt(a_y), np.sqrt(1-a_y))
            dy[j,i] = c_y * 6371229
    print(j,i)
    dy[:,i+1] = dy[:,i]
    
    return dx, dy

[dopplershift]

Are you interested in turning this into a proper PR? My only reservation about this is that it's a assuming a spherical earth. I feel like there must be a way to do this with Proj.4/pyproj.

[kgoebber]

It turns out that there is functionality in pyproj to do this very thing (see code below). I have confirmed that the previous code gives the same values when designating a spherical earth (with only changing the radius to match. I also tested other ellipsoids including WGS84 and GRS80 and the differences are at most 0.5% (one half of one percent). It appears that the NAM is based on a spherical earth to get the appropriate values at 25N and 35N.

So either way would work, it just depends on whether we would want to make a dependency on pyproj (although that might be essentially covered if we make metpy dependent on cartopy in the future).

I've got a notebook that does the calculation of distances, computes absolute vertical vorticity, and compares with the absolute vertical vorticity given directly from the model. I could make it a proper PR if we decide which path is preferable.

https://nbviewer.jupyter.org/url/fujita.valpo.edu/%7Ekgoebber/NAM_vorticity.ipynb

def new_calc_dx_dy(longitude,latitude,shape,radius=6370997.):
    ''' This definition calculates the distance 
        between grid points that are in
        a latitude/longitude format.
        
        Using pyproj GEOD; different Earth Shapes 
        https://jswhit.github.io/pyproj/pyproj.Geod-class.html
        Common shapes: 'sphere', 'WGS84', 'GRS80'
        
        Accepts, 1D arrays for latitude and longitude
        
        Returns: dx, dy; 2D arrays of distances 
                       between grid points in the x and y direction in meters 
    '''
    from pyproj import Geod
    
    if (radius != 6370997.):
        g = Geod(a=radius,b=radius)
    else:
        g = Geod(ellps=shape)
    
    dx = np.empty(latitude.shape)
    dy = np.zeros(longitude.shape)
    
    for i in range(latitude.shape[1]):
        for j in range(latitude.shape[0]-1):
            _, _, dx[j,i] = g.inv(longitude[j,i],latitude[j,i],longitude[j+1,i],latitude[j+1,i])
    dx[j+1,:] = dx[j,:]
    
    for i in range(latitude.shape[1]-1):
        for j in range(latitude.shape[0]):
            _, _, dy[j,i] = g.inv(longitude[j,i],latitude[j,i],longitude[j,i+1],latitude[j,i+1])
    dy[:,i+1] = dy[:,i]
    
    return dx, dy

[mccrayc]

I suppose any progress on this would be in part held back by numpy/numpy#9401, with numpy's gradient function requiring scalar dx and dy.

I can get a gradient field by iterating through the lat and lon coordinates of an array and calculating gradient at each individual grid point using the local dx and dy from the helper function above (see code below), but it's extremely slow (not surprising since we're performing 361*720 gradient calculations in the case of a 0.5x0.5 deg lat/lon grid. Does anyone have any ideas on a workaround for this?

field = temperature.values
gradient = np.empty([field.shape[0],field.shape[1]])
for y in np.arange(0, field.shape[0]):
    for x in np.arange(0, field.shape[1]):
        dy_point = dy[y,x]
        dx_point = dx[y,x]
        grad = np.gradient(field, dy_point.magnitude, dx_point.magnitude)            
        ddy_thta, ddx_thta = grad[-2:]
        mag_thta = np.sqrt(ddx_thta**2 + ddy_thta**2)
        gradient[y,x] = mag_thta[y,x]

[dopplershift]

I think our solution is going to be to include our own gradient function. See #174 for more info. Happy to have help there.

[dopplershift]

Fixed ini 0.7/0.8 with lat_lon_grid_spacing and now lat_lon_grid_deltas. Handling of other coordinate systems and gradient, etc. is still open in #174.