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

Question

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

Can you provide some example input data and expected output data? — Nils Werner, Nov 08 '21 at 14:47
@NilsWerner. The question is actually pretty clear about that part — Mad Physicist, Nov 08 '21 at 14:54
OP, did you copy and paste the same thing multiple times by accident? — Mad Physicist, Nov 08 '21 at 14:55
Using itertools will remove a for loop, but likely be slower. If you're using numpy arrays, vectorize — Mad Physicist, Nov 08 '21 at 14:57
Sorrry some mistake. removed duplicates now. While I was searching for answers, came across https://stackoverflow.com/questions/28684492/numpy-equivalent-of-itertools-product. As per this itertools is the fastest way to do it. But this is 4 years old. Any updated info would be of help. — arundeep78, Nov 08 '21 at 15:04

Mad Physicist · Accepted Answer · 2021-11-08T17:01:41.157

1

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)

edited Nov 08 '21 at 17:01

answered Nov 08 '21 at 15:09

Mad Physicist

107,652
25
181
264

Thanks a lot. In the meantime, I also found np.meshgrid that was much faster than itertools.product. From 4.5 sec it brought it down to 150ms. So total execution time was 350ms. i used your code ( took some time to understand :) ). It ran in my case in about 140ms. But was giving me wrong results. so I had to change it. – arundeep78 Nov 08 '21 at 16:33
I changed it to ```` grid_delta = 0.5 * grid_size # dlon: (3600,) dlon = np.full(lon.shape, grid_size / 360) # 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) return area. for whatever reason, it runs in just 43 ms !!!! It seems /360 takes rest of the time – arundeep78 Nov 08 '21 at 16:33
i accepted your answer as it answered the question about numpy martrix creation. Regarding actual result, that was not the problem statement anyways. – arundeep78 Nov 08 '21 at 16:51
I cannot edit the answer. but my units were in degrees. As per the original formula then dlon should be using np.full(lons.shape,grid_size/360) and area should be np.deg2rad(radius**2*dlat*dlon) – arundeep78 Nov 08 '21 at 17:28

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

1 Answers1