如何使用python生成常规地理网格? [英] How can I generate a regular geographic grid using python?

问题描述

我想在某个地图区域上检索规则网格的所有纬度/经度坐标对.我找到了geopy库，但根本无法解决问题.

I want to retrieve all lat/lon coordinate pairs of a regular grid over a certain map area. I have found the geopy library, but didn't manage at all to approach the problem.

例如，我有一个由纬度/经度坐标中的四个角描述的矩形地理区域，我试图计算间距为例如的网格.覆盖该区域1公里.

For example, I have a rectangular geographic area described by its four corners in lat/lon coordinates, I seek to calculate the grid with spacing of e.g. 1km covering this area.

推荐答案

初步注意事项

某些区域的定义方式略有不同.如果只是一个矩形区域(请注意:投影的矩形不一定在地球表面上是矩形的！)，则可以使用所需的步长，简单地从两个坐标维度的最小值到最大值进行迭代.如果手头有任意多边形，则需要测试哪些生成点与该多边形相交，并且仅返回满足该条件的那些坐标对.

Preliminary Considerations

It makes a slight difference how your certain area is defined. If it's just a rectangular area (note: rectangular in projection is not neccessarily rectangular on Earth's surface!), you can simply iterate from min to max in both coordinate dimensions using your desired step size. If you have an arbitrary polygon shape at hand, you need to test which of your generated points intersects with that poly and only return those coordinate pairs for which this condition holds true.

规则网格不等于投影上的规则网格.您正在谈论纬度/经度对，这是一个极坐标系，以度数表示，近似于地球表面形状.在经/纬(EPSG:4326)中，距离不是以米/公里/英里为单位，而是以度为单位.

A regular grid does not equal a regular grid across projections. You are talking about latitude/longitude pairs, which is a polar coordinate system measured in degrees on an approximation of Earth's surface shape. In lat/lon (EPSG:4326), distances are not measured in meters/kilometers/miles, but rather in degrees.

此外，我假设您要计算一个网格，其水平"步距平行于赤道(即纬度).对于其他网格(例如旋转的矩形网格，垂直于经度的垂直网格等)，您需要花费更多的精力来变换形状.

Furthermore, I assume you want to calculate a grid with its "horizontal" steps being parallel to the equator (i.e. latitudes). For other grids (for example rotated rectangular grids, verticals being parallel to longitudes, etc.) you need to spend more effort transforming your shapes.

问问自己:您要创建一个以度或米为单位的规则间隔的网格吗?

以度为单位的网格

如果要以度为单位，则可以简单地进行迭代:

If you want it in degrees, you can simply iterate:

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

但是:一定要知道，整个地球表面的长度(以米为单位)的步长(以度为单位)是不同的.例如，靠近赤道的0.001 delta纬度在表面上的距离(以米为单位)与靠近两极的距离不同.

But: Be sure to know that the length in meters of a step in degrees is not the same across Earth's surface. For example, 0.001 delta degrees of latitude close to the equator covers a different distance in meters on the surface than it does close to the poles.

以米为单位的网格

如果要以米为单位，则需要将输入区域(地图上的某些区域)的纬度/经度边界投影到支持以米为单位的距离的坐标系中.您可以使用 Haversine公式作为近似值来计算经纬度对之间的距离，但这不是您可以使用的最佳方法.

If you want it in meters, you need to project your lat/lon boundaries of your input area (your certain area on a map) into a coordinate system that supports distances in meters. You can use the Haversine formula as a rough approximation to calculate the distance between lat/lon pairs, but this is not the best method you can use.

更好的方法是搜索合适的投影，将您感兴趣的区域转换为该投影，通过简单的迭代创建网格，获取点，然后将其投影回经纬度对.例如，欧洲的合适投影将是EPSG:3035.顺便说一下，Google Maps将EPSG:900913用于其网络地图服务.

Better yet is to search for a suitable projection, transform your area of interest into that projection, create a grid by straightforward iteration, get the points, and project them back to lat/lon pairs. For example, a suitable projection for Europe would be EPSG:3035. Google Maps uses EPSG:900913 for their web mapping service, by the way.

在python中，您可以使用库shapely和pyproj处理地理形状和投影:

In python, you can use the libraries shapely and pyproj to work on geographic shapes and projections:

import shapely.geometry
import pyproj

# Set up projections
p_ll = pyproj.Proj(init='epsg:4326')
p_mt = pyproj.Proj(init='epsg:3857') # metric; same as EPSG:900913

# 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 = pyproj.transform(p_ll, p_mt, sw.x, sw.y) # Transform NW point to 3857
transformed_ne = pyproj.transform(p_ll, p_mt, 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(pyproj.transform(p_mt, p_ll, 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))

此示例生成此均匀间隔的网格:

This example generates this evenly-spaced grid:

这篇关于如何使用python生成常规地理网格?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！