如何计算每个lat/long网格正方形的m2大小

Question

我有数据的单位m2的网格分辨率的估计。我需要计算每个纬度/经度网格单元格中的m2数？

在两极附近，细胞尺寸比赤道小得多，所以这一点很重要。

我想要一个netcdf文件或数组，在每个网格平方中的平方米数。

Answer 1

如果有人想要一个在每个拉长网格单元中的平方公尺数的网卡。这可能不是最干净的解决方案，但将使用xarray在每个网格中创建一个m2的netcdf (M2)。

gridsize()函数来自另一个堆栈溢出问题。然后，我们可以制作一个虚拟数组，并使用每个位置的经度距离创建一个m2s地球场。

"""
This will create a global grid of the approximate size of each grid square.
"""
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

def gridsize(lat1):
   #https://en.wikipedia.org/wiki/Haversine_formula
   #https://stackoverflow.com/questions/639695/how-to-convert-latitude-or-longitude-to-meters/11172685#11172685
   lon1=200
   import math
   lat2=lat1
   lon2=lon1+1

   R = 6378.137 # // Radius of earth in km
   dLat = lat2 * np.pi / 180 - lat1 * np.pi / 180
   dLon = lon2 * np.pi / 180 - lon1 * np.pi / 180
   a = np.sin(dLat/2) * np.sin(dLat/2) + np.cos(lat1 * np.pi / 180) * np.cos(lat2 * np.pi / 180) * np.sin(dLon/2) * np.sin(dLon/2)
   c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))
   d = R * c
   return d * 1000 #; // meters


boxlo,boxla=np.array(np.meshgrid(np.arange(-179.5,179.5,1),np.arange(-89.5,89.5,1)))
sizes=np.ones(boxlo.shape)
grid=gridsize(boxla)

grid_nc=xr.DataArray(grid,coords={'lat':boxla[:,1],'lon':boxlo[1,:]},dims=['lat','lon'])
lat_size=110567 #in m
grid_nc['m2']=grid_nc*lat_size
grid_nc=grid_nc['m2']
grid_nc.to_netcdf('earth_m2.nc')

plt.pcolormesh(boxlo[1,:],boxla[:,1],grid_nc)
plt.colorbar()
plt.show()

Answer 2

用球面三角学方法，半径R的球面上的三角形在北极有1个顶点(纬度:π/2)，在同一纬度上有2个顶点- -π/2

S(x, d) = (cos⁻¹((cos(b) - cos²(a))/sin²(a)) + 2cos⁻¹((cos(a) - cos(a)cos(b))/(sin(a)sin(b))) - π)R² where a = R(π/2 - x) and b = Rd

因此，半径为R的球面上的网格矩形在d半径与纬度x₁> x₂相隔的经线之间的表面积是

S(x₂, d) - S(x₁, d)

Answer 3

一种选择是将你的细胞转换成一个坐标参考系统(CRS)，它的单位是米，而不是度。这样，面积的计算就简单了。

我猜你的坐标在WGS84。

对于目标CRS，有一些选择，特别是当您知道点的位置时，但是像这样的全局CRSs的一个公共集合是通用横向墨卡托(UTM)，或者靠近通用极地平片极点。

例如，对于UTM，假设表单[lon, lat]的点列表中的最后一点等于第一个点：

import pyproj
from shapely.geometry import Polygon
from shapely.ops import transform

def utm_epsg(lon: float, lat: float) -> int:
  """
  Return the UTM EPSG code for the given lon-lat.
  """
  offset = int(round((183 + lon) / 6.0))
  return 32600 + offset if lat > 0 else 32700 + offset

for lat in range(-79, 83):
    for lon in range(-179, 179):
        polygon = Polygon([
            [lon, lat],
            [lon+1, lat],
            [lon+1, lat+1],
            [lon, lat+1],
            [lon, lat],
        ])
        
        src_crs = pyproj.CRS.from_epsg(4326)
        tgt_crs = pyproj.CRS.from_epsg(utm_epsg(polygon.centroid.x, polygon.centroid.y))

        project = pyproj.Transformer.from_crs(src_crs, tgt_crs, always_xy=True).transform

        utm_polygon = transform(project, polygon)

        # aggregate into some result.  Here just printed to stdout.
        print(polygon.centroid, utm_polygon.area)

值得注意的是，UTM的定义不是80°S以南，84°N以北。