将Shapefile转换为2D网格稀疏矩阵

Question

我在使用NumPy进行地理空间计算时，被所有附加库搞得眼花缭乱。

如何以最简单的方式获取一个覆盖整个地球范围的shapefile，并从中构建一个稀疏矩阵，该矩阵表示具有指定分辨率的纬度-经度网格，在shapefile多边形内的所有网格点矩阵条目为1，其他地方为0？

我知道如何使用GeoPandas或Fiona读取shapefile；但是我卡在了“转换为稀疏矩阵”的部分。看起来“栅格化”需要另一个附加库，无法控制网格分辨率，并且无法输出除TIFF之外的任何格式，这不太有用。

- zwol

---

你需要多次这样做，还是可以使用GIS只做一次并将该栅格输出导入Python / NumPy？ - Ian Turton

我只需要做一次，但是到目前为止，我使用独立的GIS的经验是“它们都不起作用”，因此我更倾向于使用NumPy和相关工具。 - zwol

你可以在Gis.stackexchange.com上提问。 - Ian Turton

Answer 1

这里有一种方法：

1. 创建一个空的稀疏矩阵
2. 编写一个函数将随机点发送到“小”点集
3. 对于每个多边形，在包围盒中均匀采样所有点。对于每个落入多边形的点，应用上述函数以获得标准的key，并使用1的键和值更新稀疏矩阵。

我不确定它是否属于简单类别，但可以使用osgeo和scipy在python中完成。当然，如果您有大型多边形，则采样速度将非常慢，但由于您使用了稀疏矩阵，因此我认为这不会成为问题。您也可以在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()

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

- hilberts_drinking_problem