我在使用NumPy进行地理空间计算时,被所有附加库搞得眼花缭乱。
如何以最简单的方式获取一个覆盖整个地球范围的shapefile,并从中构建一个稀疏矩阵,该矩阵表示具有指定分辨率的纬度-经度网格,在shapefile多边形内的所有网格点矩阵条目为1,其他地方为0?
我知道如何使用GeoPandas或Fiona读取shapefile;但是我卡在了“转换为稀疏矩阵”的部分。看起来“栅格化”需要另一个附加库,无法控制网格分辨率,并且无法输出除TIFF之外的任何格式,这不太有用。
我在使用NumPy进行地理空间计算时,被所有附加库搞得眼花缭乱。
如何以最简单的方式获取一个覆盖整个地球范围的shapefile,并从中构建一个稀疏矩阵,该矩阵表示具有指定分辨率的纬度-经度网格,在shapefile多边形内的所有网格点矩阵条目为1,其他地方为0?
我知道如何使用GeoPandas或Fiona读取shapefile;但是我卡在了“转换为稀疏矩阵”的部分。看起来“栅格化”需要另一个附加库,无法控制网格分辨率,并且无法输出除TIFF之外的任何格式,这不太有用。
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()
这是一张图片;我省略了一些大国家以加快速度。