Python中地理坐标的区域元素

0 人关注

我需要对网格化的天气数据（降雨量）进行一些整合/汇总，以标准的地理坐标。因此，我需要计算与我的每个网格点相关的面积元素。基本上是网格单元的实际面积。我们谈论的是一个标准的、等距的、直线型的纬度/纬度网格。

我怎样才能在Python中做到这一点？

我知道我可以用球面坐标近似法合理地获得它们，但这必须是超级标准的，所以不妨使用一个有适当地球形状模型的软件包。然而，我找了很多地方，都没有发现如何在普通的地理空间包中完成这一基本任务。似乎它通常隐藏在引擎盖或其他地方，在内部处理。但我相信可以通过某种方式把它挖出来。

编辑：这里有一个手工制作的解决方案。如何计算每个经纬度方格的面积，以平方米计，但如果可能的话，我宁愿使用标准库。

python

integration

geospatial

coordinate-systems

coordinate-transformation

Ben Farmer

发布于 2022-06-16

2 个回答

Olivier Archer

发布于 2022-06-17

已采纳

0 人赞同

See pyproj.Geod.inv ，可用于返回两点之间的距离（以米为单位）（以lon/lat为单位）。

g=Geod(ellps='WGS84')
lon2D,lat2D = np.meshgrid(np.arange(0,20,0.1),np.arange(30,45,0.1))
_,_, distEW = g.inv(lon2D[:,:-1],lat2D[:,1:], lon2D[:,1:], lat2D[:,1:])
_,_, distNS = g.inv(lon2D[1:,:],lat2D[1:,:], lon2D[1:,:], lat2D[:-1,:])
pixel_area = distEW[1:,:] * distNS[:,1:]


         0
         
         人赞同


          
           奥利维尔的答案很好，也很简洁，但为了我自己的好奇心，我想看看它与更精确的计算相比如何（他们的答案假设网格单元足够小，可以近似为正方形）。
          
          import numpy as np
from pyproj import Geod
from shapely.geometry import LineString, Point, Polygon
lons = np.arange(0,20,0.1)
lats = np.arange(30,45,0.1)
# Fast way, square approx of grid cells
geod = Geod(ellps="WGS84")
lon2D,lat2D = np.meshgrid(lons, lats)
_,_, distEW = geod.inv(lon2D[:,:-1],lat2D[:,1:], lon2D[:,1:], lat2D[:,1:])
_,_, distNS = geod.inv(lon2D[1:,:],lat2D[1:,:], lon2D[1:,:], lat2D[:-1,:])
square_area = distEW[1:,:] * distNS[:,1:]
# Slower way, geodesic areas
def stack_coords(x, y):
    """Stacks coordinate lists into a 2D array of coordinate pairs,
       flattened to list of coordinate pairs"""
    out = np.vstack(np.stack(np.meshgrid(x, y), axis=2))
    return out
def get_lat_squares(lats, lons):
    """Construct shapely Polygons for a single longitude but
       every latitude.
       Area is constant with longitude so just copy results.
    lats_0 = lats[:-1]
    lats_1 = lats[1:]
    lons_0 = lons[0:1]
    lons_1 = lons[1:2]
    c1 = stack_coords(lons_0, lats_0)
    c2 = stack_coords(lons_1, lats_0)
    c3 = stack_coords(lons_1, lats_1)
    c4 = stack_coords(lons_0, lats_1)
    squares = []
    for p1, p2, p3, p4 in zip(c1, c2, c3, c4):
        squares.append(Polygon([p1, p2, p3, p4]))
    return squares
def get_areas(lats, lons):
    squares = get_lat_squares(lats, lons)
    geod = Geod(ellps="WGS84")
    areas = []
    for square in squares:
        area = geod.geometry_area_perimeter(square)[0]
        areas.append(area)
    return np.array(areas)
geodesic_areas = get_areas(lats, lons)
for a1, a2 in zip(geodesic_areas, square_area[:,0]):
   print(a1, a2)
Output:
106904546.2618103 106850809.52460858
106798611.31295013 106744711.14938568
106692351.02400208 106638287.58503169
106585765.66596985 106531539.10307406
106478855.51091766 106424465.9760592
88300037.70746613 88224423.89253764