Shapefile zu 2D-Gitter als spärliche Matrix

Question

Ich bin in all den Bolt-on-Bibliotheken für die Durchführung von Geodatenberechnungen mit NumPy verloren gegangen.

Was ist der direkteste Weise eine Formdatei zu nehmen (dessen Umfang ist die gesamte Erde), und aus ihr eine Sparse-Matrix-Konstrukt, das ein Breiten-Längen-Raster mit einer spezifizierten Auflösung darstellt, mit der Matrix Einträge 1 überhaupt Rasterpunkte innerhalb der Polygone der Formdatei, 0 woanders?

Ich kann in einem Shapefile mit GeoPandas oder Fiona lesen; Wo ich feststecke, ist der Teil, der in eine dünne Matrix umgewandelt wird. "Rasterisierung" scheint noch einen weiteren Anstoß zu erfordern, gibt mir nicht die Kontrolle über die Rasterauflösung und kann nichts anderes als TIFFs ausspucken, was weniger als nützlich ist. Hier

Answer 1

ist ein Ansatz:

1. Erstellen leer Sparse Matrix
2. eine Funktion schreiben beliebigen Punkt auf einer „kleinen“ Menge von Punkten
3. Für jedes Polygon zu schicken, probieren alle Punkte in der Bounding Box einheitlich. Für jeden Punkt, der in die Polygon fällt, gelten obige Funktion Standard key zu erhalten und die spärliche Matrix mit dem Schlüssel und den Wert der 1.

ich nicht sicher aktualisieren, wenn es unter einfacher Kategorie fällt, aber es kann getan werden, in Python mit osgeo und scipy. Natürlich ist die Abtastung sehr langsam, wenn Sie große Polygone haben, aber da Sie eine dünne Matrix verwenden, nehme ich an, dass dies kein Problem ist. Sie können die Auflösung anpassen und mit Projektionen innerhalb von osgeo spielen.

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() 

Hier ist ein Bild; Ich habe einige große Länder ausgelassen, um die Dinge zu beschleunigen.