我有固定解析度的經度和緯度陣列，即 .1。這給了我 1800 lats 和 3600 lons。我想創建一個 1800x 3600 的矩陣，它將根據此處的公式存盤每個網格的區域。IE
A = 2 pi R^2 |sin(lat1)-sin(lat2)| |lon1-lon2|/360
我有 lons 已經在代表網格中心的陣列中。目前我使用一個公式，它計算給定矩形框的面積。
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
我在每個緯度/經度組合的雙回圈中使用它來獲取 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)
正如預期的那樣，這很慢。我不確定我可以使用哪種方法來提高效率。我瀏覽了論壇以創建 NxM 矩陣，但無法獲得此問題的解決方案。
在寫這個問題時，在 stackoverflow 上遇到了這個執行緒來使用 itertools.chain。如果有幫助，將嘗試按照此更改我的代碼。將更新我的發現。
與此同時，任何朝著正確方向的幫助都會有所幫助。
更新：我使用 itertools.product 更改了我的代碼
lat_longs = np.array(list(itertools.product(*[lats.tolist(),longs.tolist()])))
并更新函式以接受質心。
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
然后我使用重新排列回傳區域陣列
areas_mat = areas.reshape((lats.shape[0], longs.shape[0]))
現在代碼中最長的部分是 itertools.product。大約需要4.5秒，而面積計算只需要大約350ms。
還有其他方法可以更快地獲得第一個組合嗎？
uj5u.com熱心網友回復：
您可以在所有陣列中輕松矢量化此操作。給定一個latsshape(1800,)的陣列和一個shape的陣列，您可以對它們lons進行整形(3600,)，以便廣播計算產生一個正確形狀的陣列。
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)
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/352944.html
上一篇：在imshow上疊加等高線圖
下一篇：定期切片2Dnumpy陣列