Odpowiedzi:
Python GDAL / OGR Cookbook ma jakiś przykładowy kod do buforu geometrię .
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())
oraz do obliczenia przecięcia między dwiema geometriami
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()
Geometria może być odczytywana i zapisywana w plikach kształtu i wielu innych formatach.
Aby uprościć, Shapely: manual pozwala na całą obróbkę geometryczną PostGIS w Pythonie.
Pierwszą przesłanką Shapely jest to, że programiści Python powinni mieć możliwość wykonywania operacji geometrii typu PostGIS poza RDBMS ...
Pierwszy przykład 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), ....}
Przykład wielokąta z 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) -> GeoJSON
Drugą przesłanką jest to, że trwałość, serializacja i odwzorowanie cech są znaczącymi, ale ortogonalnymi problemami. Możesz nie potrzebować stu czytników i pisarzy formatu GIS lub mnóstwa projekcji State Plane, a Shapely ich nie obciąży.
Więc łączysz go z innymi modułami Python, aby czytać lub zapisywać pliki kształtu i manipulować projekcjami jak osgeo.ogr, Fiona lub PyShp .
Przeszukując Gis StackExchange, możesz znaleźć wiele przykładów, ale podam jeszcze jeden, aby zilustrować połączenie foremnej i Fiony oraz użycie funkcji forex przecięcia () i buffer () (Można to zrobić za pomocą PyShp).
Biorąc pod uwagę dwa pliki kształtów polilinii:
Oblicz przecięcie (funkcja przecięcia () z foremnego)
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)
Zapisz wynik jako nowy plik kształtu
# 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}})
Wynik:
Buforuj pojedyncze punkty (bufor funkcji () foremnego)
# 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}})
Wynik
Buforuj geometrię MultiPoint
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}})
Moją biblioteką geo-przetwarzania „przejdź do” jest „Biblioteka teledetekcji i GIS” (RSGISLib). Jest łatwy w instalacji i obsłudze, a dokumentacja jest naprawdę dobra. Ma funkcjonalność przetwarzania wektorowego i rastrowego - bardzo rzadko muszę zbliżać się do GUI. Można go znaleźć tutaj: http://rsgislib.org .
Przykładem w tym przypadku jest:
rsgislib.vectorutils.buffervector(inputvector, outputvector, bufferDist, force)
Polecenie buforowania wektora o określoną odległość.
Gdzie:
Przykład:
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)