Create NxM array from 2 arrays/lists of N and M size-CodePudding

I have longitude and latitude arrays of fixed resolutions i.e. .1. This gives me 1800 lats and 3600 lons. I want to create a matrix of 1800x 3600 that will store area for each grid based on the formula here . i.e.

A = 2piR^2 |sin(lat1)-sin(lat2)| |lon1-lon2|/360

I have lons are lats already in arrays which represents centre of the grid. Currently I use a formula, which calculates area for a given rectangle box.

def grid_area(lat1, lon1, lat2, lon2, radius= 6365000):
    """
    Calculate grid area based on lat-long points of rectangle/square grid size by degrees.
    Calculations are without any prohection system.
    
    radius in meters is used to make it generic. Defaults to Earth 
    
    Formuala from : https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2004/msg00023.html
    """
    import numpy as np
    area = (np.pi/180)*(radius**2) *np.abs(np.sin(np.radians(lat1)) - np.sin(np.radians(lat2))) * np.abs(lon1 -lon2)/360
    return area

I use this in a double loop for each lat/lon combination to get the area_grid.

  grid_areas = np.zeros((len(lats), len(longs)))
  for ll in range(len(longs)-1):
     for lt in range(len(lats)-1):
        lt1 = np.round(lats[lt] .05,2)
        ll1 = np.round(longs[ll]-.05,2)
        lt2 = np.round(lats[lt]-.05,2)
        ll2 = np.round(longs[ll] .05,2)
        
        grid_areas[lt,ll] = grid_area(lt1,ll1,lt2,ll2)

This as expected is slow. I am not sure which approach I can use to make it efficient. I looked through the forum to create NxM matrixes, but not able to get the solution for this problem.

While writing this question, came across this thread on stackoverflow to use itertools.chain. Will try to change my code as per this, if that helps. Will update my findings on that.

In the meantime, any help in the right direction would help.

UPDATE: I changed my code using itertools.product

lat_longs = np.array(list(itertools.product(*[lats.tolist(),longs.tolist()])))

and updated the function to accept centroids.

def grid_area(lat=None, lon=None, grid_size=.1, radius= 6365000):
    """
    Calculate grid area based on lat-long points of rectangle/square grid size by degrees.
    Calculations are without any prohection system.
    
    radius in meters is used to make it generic. Defaults to Earth 
    
    Formuala from : https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2004/msg00023.html
    """
    import numpy as np
    grid_delta = grid_size/2
    lat1 = lat grid_delta
    lat2 = lat-grid_delta
    
    lon1 = lon - grid_delta
    lon2 = lon   grid_delta
    
    area = (np.pi/180)*(radius**2) *np.abs(np.sin(np.radians(lat1)) - np.sin(np.radians(lat2))) * np.abs(lon1 -lon2)/360
    return area

I then rearrange the return area array using

areas_mat = areas.reshape((lats.shape[0], longs.shape[0]))

Now the longest part of the code is the itertools.product. it takes about 4.5 seconds, while the area calculation takes only about 350ms.

Any other way to get that first combination faster?

Update2: Final code

Once I tried, I found that area was not correct, even when the code was aligned with formula in the link. used the 2nd source for final version. Final code is

def grid_area_vec(lat=None, lon=None, grid_size=.1, radius= 6365000):
    """
    Calculate grid area based on lat-long points of rectangle/square grid size by degrees.
    Calculations are without any prohection system.
    
    radius in meters is used to make it generic. Defaults to Earth 
    
    Orig Formula from : https://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2004/msg00023.html
    Another source for formula, finally used
    https://gis.stackexchange.com/questions/413349/calculating-area-of-lat-lon-polygons-without-transformation-using-geopandas
    
    """
    import numpy as np
    
    grid_delta = 0.5 * grid_size
    # dlon: (3600,)
    dlon = np.full(lon.shape, np.deg2rad(grid_size))
    # dlat: (1800, 1)
    dlat = np.abs(np.sin(np.deg2rad(lat   grid_delta)) -
                  np.sin(np.deg2rad(lat - grid_delta)))[:, None]
    # area: (1800, 3600)
#     area = np.deg2rad(radius**2 * dlat * dlon)
    area = radius**2 * (dlat * dlon)
    
    return area

CodePudding user response：

You can trivially vectorize this operation across all your arrays. Given an array lats with shape (1800,), and an array lons with shape (3600,), you can reshape them so that the broadcasted computation yields an array of the correct shape.

grid_delta = 0.5 * grid_size
# dlon: (3600,)
dlon = np.full(lons.shape, np.rad2deg(grid_size))
# dlat: (1800, 1)
dlat = np.abs(np.sin(np.deg2rad(lats   grid_delta)) -
              np.sin(np.deg2rad(lats - grid_delta)))[:, None]
# area: (1800, 3600)
area = np.rad2deg(radius**2 * dlat * dlon)