如何使用Python生成规则的地理网格？

Question

我想获取特定地图区域上的所有经纬度坐标对，以形成一个常规网格。我已经找到了geopy库，但是没有成功解决这个问题。

例如，我有一个由其四个角落的经纬度坐标描述的矩形地理区域，我希望计算出覆盖该区域的间距为1公里的网格。

- Amira Akhdal

---

类似于NSEW的基本点？如何在(lon, lat)对中给出基本点，并且它们如何具有距离？请详细阐述您的问题并添加到目前为止编写的代码。 - jbndlr

抱歉，我错了。我想要找到某个区域内的所有经纬度坐标。例如，在谷歌地图上，我有一个区域的4个点的经纬度，我想检索该区域内所有相距10公里的经纬度坐标。我现在还没有代码，我仍在使用geopy进行搜索。 - Amira Akhdal

你想要在给定的多边形内创建一个点坐标的正则网格吗？ - jbndlr

是的，这就是！在geopy或其他库中有什么东西吗？ - Amira Akhdal

Answer 1

初步考虑

如何定义您的特定区域略有不同。如果它只是一个矩形区域（注意：在投影中是矩形的，但在地球表面上不一定是），您可以在两个坐标维度中使用所需的步长从最小值到最大值进行迭代。如果您手头有任意多边形形状，则需要测试哪些生成的点与该多边形相交，并仅返回满足此条件的坐标对。

计算规则网格

在不同的投影中，规则网格并不相等。您正在谈论纬度/经度对，这是以度数为单位测量地球表面形状近似的极坐标系。在lat/lon（EPSG:4326）中，距离不是用米/千米/英里来衡量的，而是用度数来衡量。

此外，我假设您希望计算其“水平”步长与赤道（即纬度）平行的网格。对于其他网格（例如旋转矩形网格，垂直线与经度平行等），您需要花费更多的精力来转换形状。

问问自己：您想要在度数还是米中创建等间距网格？

以度为单位的网格

如果您想要以度为单位的网格，可以简单地迭代：

stepsize = 0.001
for x in range(lonmin, lonmax, stepsize):
    for y in range(latmin, latmax, stepsize):
        yield (x, y)

但是: 请注意，一度的步长长度在地球表面上的米数并不相同。例如，接近赤道的0.001度纬度与接近极地的0.001度纬度所覆盖的距离在地图上显示的长度是不同的。

米制网格

如果你需要使用米制单位，你需要将你所关注的特定区域（地图上的）的纬度/经度范围投影到支持使用米制单位测量距离的坐标系中。你可以使用Haversine公式来粗略地计算纬度/经度对之间的距离，但这并不是最好的方法。

更好的做法是寻找一个合适的投影方式，将你感兴趣的区域转换为该投影方式，通过简单迭代创建一个网格，获取点，并将它们重新投影回纬度/经度对。例如，欧洲的一个适合的投影方式是EPSG：3035。顺便说一下，Google地图使用EPSG：900913作为他们的Web地图服务的投影方式。

在Python中，您可以使用shapely和pyproj库来处理地理形状和投影：

import shapely.geometry
import pyproj

# Set up transformers, EPSG:3857 is metric, same as EPSG:900913
to_proxy_transformer = pyproj.Transformer.from_crs('epsg:4326', 'epsg:3857')
to_original_transformer = pyproj.Transformer.from_crs('epsg:4326', 'epsg:3857')

# Create corners of rectangle to be transformed to a grid
sw = shapely.geometry.Point((-5.0, 40.0))
ne = shapely.geometry.Point((-4.0, 41.0))

stepsize = 5000 # 5 km grid step size

# Project corners to target projection
transformed_sw = to_proxy_transformer.transform(sw.x, sw.y) # Transform NW point to 3857
transformed_ne = to_proxy_transformer.transform(ne.x, ne.y) # .. same for SE

# Iterate over 2D area
gridpoints = []
x = transformed_sw[0]
while x < transformed_ne[0]:
    y = transformed_sw[1]
    while y < transformed_ne[1]:
        p = shapely.geometry.Point(to_original_transformer.transform(x, y))
        gridpoints.append(p)
        y += stepsize
    x += stepsize

with open('testout.csv', 'wb') as of:
    of.write('lon;lat\n')
    for p in gridpoints:
        of.write('{:f};{:f}\n'.format(p.x, p.y))

这个例子生成了如下等间隔的网格:

- jbndlr

---

非常感谢，您讲解得非常清晰，教给了我很多新的东西，谢谢！ - Amira Akhdal

你的代码有一个小错误：在写入 testout.csv 时，流被定义为“wb”，但是却写入了一个字符串。应该将其更改为“w”。 - Aleksandar Jovanovic

你混淆了角落的名称。nw 不是左上角，而是左下角，因此应该是 sw。同样，se 应该被命名为 ne，表示东北角（右上角）。 - dasjanik

@Badmiral 我在 QGIS 中加载了生成的 CSV 并截取了屏幕截图，这对我来说是最快的解决方案。你也可以轻松使用 matplotlib 等工具。 - jbndlr

很好的解决方案，不过有一个小错误。应该修改为：to_original_transformer = pyproj.Transformer.from_crs('epsg:3857', 'epsg:4326') - BSnider