Gibt es mit Ausnahme von ArcPy eine Python-Bibliothek, die Geoverarbeitung mit Shapefiles durchführen kann, z. B. Buffer / Intersect?
Das Python GDAL / OGR-Kochbuch enthält einen Beispielcode zum Puffern einer Geometrie .
from osgeo import ogr
wkt = "POINT (1198054.34 648493.09)"
pt = ogr.CreateGeometryFromWkt(wkt)
bufferDistance = 500
poly = pt.Buffer(bufferDistance)
print "%s buffered by %d is %s" % (pt.ExportToWkt(), bufferDistance, poly.ExportToWkt())und zum Berechnen des Schnittpunkts zwischen zwei Geometrien
from osgeo import ogr
wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"
poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)
intersection = poly1.Intersection(poly2)
print intersection.ExportToWkt()Die Geometrien können in Shapefiles und einer Vielzahl anderer Formate gelesen und geschrieben werden .
Zur Vereinfachung erlaubt Shapely: manual die gesamte Geometrieverarbeitung von PostGIS in Python.
Die erste Prämisse von Shapely ist, dass Python-Programmierer in der Lage sein sollten, Geometrieoperationen vom Typ PostGIS außerhalb eines RDBMS auszuführen ...
Das erste Beispiel für PolyGeo
from shapely.geometry import Point, LineString, Polygon, mapping
from shapely.wkt import loads
pt = Point(1198054.34,648493.09)
# or
pt = loads("POINT (1198054.34 648493.09)")
bufferDistance = 500
poly = pt.buffer(bufferDistance)
print poly.wkt
'POLYGON ((1198554.3400000001000000 648493.0899999999700000, 1198551.9323633362000000
# GeoJSON
print mapping(poly)
{'type': 'Polygon', 'coordinates': (((1198554.34, 648493.09), (1198551.9323633362, 648444.0814298352), (1198544.7326402017, 648395.544838992), ....}Das Beispiel des Polygons von PolyGeo:
poly1 = Polygon([(1208064.271243039,624154.6783778917), (1208064.271243039,601260.9785661874), (1231345.9998651114,601260.9785661874),(1231345.9998651114,624154.6783778917),(1208064.271243039,624154.6783778917)])
poly2 = loads("POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"
intersection = poly1.intersection(poly2)
print intersection.wkt
print mapping(intersection) -> GeoJSONDie zweite Voraussetzung ist, dass die Persistenz, Serialisierung und Kartenprojektion von Merkmalen signifikante, aber orthogonale Probleme sind. Möglicherweise benötigen Sie keine hundert GIS-Format-Leser und -Schreiber oder die Vielzahl von State-Plane-Projektionen, und Shapely belastet Sie nicht damit.
Sie können es also mit anderen Python-Modulen kombinieren, um Shapefiles zu lesen oder zu schreiben und Projektionen wie osgeo.ogr, Fiona oder PyShp zu bearbeiten .
Wenn Sie in Gis StackExchange suchen, finden Sie viele Beispiele, aber ich gebe Ihnen ein anderes, um die Kombination von shapely und Fiona und die Verwendung der shapely-Funktionen intersection () und buffer () zu veranschaulichen (Dies hätte mit PyShp geschehen können).
Gegeben zwei Polylinien-Shapefiles:
Berechnen Sie die Schnittmenge (Funktion intersection () von shapely)
from shapely.geometry import Point, Polygon, MultiPolygon, MumtiPoint, MultiLineString,shape, mapping
import fiona
# read the shapefiles and transform to MultilineString shapely geometry (shape())
layer1 = MultiLineString([shape(line['geometry']) for line in fiona.open('polyline1.shp')])
layer2 = MultiLineString([shape(line['geometry']) for line in fiona.open('polyline2.shp')])
points_intersect = layer1.intersection(layer2)Speichern Sie das Ergebnis als neues Shapefile
# schema of the new shapefile
schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
# write the new shapefile (function mapping() of shapely)
with fiona.open('intersect.shp','w','ESRI Shapefile', schema) as e:
e.write({'geometry':mapping(points_intersect), 'properties':{'test':1}})Ergebnis:
Einzelne Punkte puffern (Funktionspuffer () von Shapely)
# new schema
schema = {'geometry': 'Polygon','properties': {'test': 'int'}}
with fiona.open('buffer.shp','w','ESRI Shapefile', schema) as e:
for point in points:
e.write({'geometry':mapping(point.buffer(300)), 'properties':{'test':1}})Ergebnis
Puffern Sie die MultiPoint-Geometrie
schema = {'geometry': 'MultiPolygon','properties': {'test': 'int'}}
points.buffer(300)
with fiona.open('buffer2.shp','w','ESRI Shapefile', schema) as e:
e.write({'geometry':mapping(points.buffer(300)), 'properties':{'test':1}})
Formschön ermöglicht Python den Zugriff auf GEOS, das Puffer / Überschneidungen / etc. Ausführen kann. GEOS ist die Bibliothek, mit der die meisten OSGeo-Programme diese Vorgänge ausführen.
Hier ist meine Liste der Python-Geoverarbeitungssoftware.
Meine Geoverarbeitungsbibliothek ist die Fernerkundungs- und GIS-Bibliothek (RSGISLib). Es ist einfach zu installieren und zu benutzen und die Dokumentation ist wirklich gut. Es bietet Funktionen für die Vektor- und Rasterverarbeitung - ich muss mich nur sehr selten einer GUI nähern. Es kann hier gefunden werden: http://rsgislib.org .
Ein Beispiel in diesem Fall ist:
rsgislib.vectorutils.buffervector(inputvector, outputvector, bufferDist, force)Ein Befehl zum Puffern eines Vektors um einen bestimmten Abstand.
Wo:
Beispiel:
from rsgislib import vectorutils
inputVector = './Vectors/injune_p142_stem_locations.shp'
outputVector = './TestOutputs/injune_p142_stem_locations_1mbuffer.shp'
bufferDist = 1
vectorutils.buffervector(inputVector, outputVector, bufferDist, True)