Shapefile 到 2D 网格作为稀疏矩阵答案

【问题标题】：Shapefile to 2D grid as sparse matrixShapefile 到 2D 网格作为稀疏矩阵
【发布时间】：2016-03-30 19:37:04
【问题描述】：

我已经完全沉浸在使用 NumPy 进行地理空间计算的所有附加库中了。

获取 shapefile（其范围是整个地球）并从中构造表示具有指定分辨率的纬度-经度网格的稀疏矩阵的最直接方法是什么？ shapefile 多边形内所有网格点的条目 1，其他地方为 0？

我知道如何使用 GeoPandas 或 Fiona 读取 shapefile；我卡住的地方是“转换为稀疏矩阵”部分。 “光栅化”似乎需要另一个附加功能，不能让我控制网格分辨率，并且除了 TIFF 之外不能吐出任何东西，这没什么用。

【问题讨论】：

您是否需要多次执行，或者您可以使用 GIS 执行一次，然后将该栅格输出导入 python/numpy？
我只需要做一次，但到目前为止我对独立 GIS 的体验是“它们都不起作用”，所以我更喜欢坚持使用 NumPy 和朋友。
你可以在 Gis.stackexchange.com 上询问

标签： python numpy geospatial sparse-matrix shapefile

【解决方案1】：

这是一种方法：

创建空稀疏矩阵
编写一个函数，将一个随机点发送到一组“小”点
对于每个多边形，统一采样边界框中的所有点。对于落入 Polygon 的每个点，应用上述函数以获取标准 key 并使用键和值 1 更新稀疏矩阵。

我不确定它是否属于简单的类别，但可以在 python 中使用osgeo 和scipy 完成。当然，如果您有大多边形，采样会非常慢，但是由于您使用的是稀疏矩阵，我认为这不会成为问题。您也可以在 osgeo 中调整分辨率和玩投影。

from itertools import product
from scipy.sparse import dok_matrix
import numpy as np

# https://pcjericks.github.io/py-gdalogr-cookbook
from osgeo import ogr

# DATA:
# http://www.naturalearthdata.com/downloads/110m-cultural-vectors/

SHP_FNAME = 'ne_110m_admin_0_countries.shp'

driver = ogr.GetDriverByName('ESRI Shapefile')
data = driver.Open(SHP_FNAME, 0)
layer = data.GetLayer()

XDIAM = 360.0
YDIAM = 180.0
XRES = YRES = 10 ** 2
dX = XDIAM / XRES
dY = YDIAM / YRES

def to_key(pt):
    x, y = pt
    x -= x % dX - XDIAM / 2
    y -= y % dY - YDIAM / 2
    return (x / dX, y / dY)

def geom_to_keys(g):
    xmin, xmax, ymin, ymax = g.GetEnvelope()
    print xmax, ymax, xmin, ymin
    xs = np.linspace(xmin, xmax, (xmax - xmin) / dX)
    ys = np.linspace(ymin, ymax, (ymax - ymin) / dY)
    for x, y in product(xs, ys):
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(x, y)
        if g.Contains(point):
            yield to_key((x, y))

smatrix = dok_matrix((XRES + 1, YRES + 1), np.int8)

one = np.int8(1)

for feature in layer:
    geom = feature.GetGeometryRef()
    if geom.Area() > 1000:
        continue
        # sampling is slow for large polygons

    for key in geom_to_keys(geom):
        smatrix.update({
            key : one,
            })

if XRES * YRES < 10 ** 6 + 1:
    from matplotlib import pyplot as plt
    plt.pcolor(smatrix.toarray().transpose())
    plt.show()

这是一张图片；为了加快速度，我省略了一些大国。

【讨论】：