Excluyendo ArcPy, ¿hay alguna biblioteca de Python que pueda hacer geoprocesamiento, como buffer / intersect, con shapefiles?
El libro de cocina Python GDAL / OGR tiene un código de muestra para almacenar una geometría .
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())y para calcular la intersección entre dos geometrías
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()Las geometrías se pueden leer y escribir en archivos de forma y en una variedad de otros formatos.
Para simplificar, Shapely: manual permite todo el procesamiento de geometría de PostGIS en Python.
La primera premisa de Shapely es que los programadores de Python deberían poder realizar operaciones de geometría de tipo PostGIS fuera de un RDBMS ...
El primer ejemplo de 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), ....}El ejemplo del polígono de 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) -> GeoJSONLa segunda premisa es que la persistencia, la serialización y la proyección del mapa de características son problemas significativos, pero ortogonales. Es posible que no necesite cientos de lectores y escritores en formato GIS o la multitud de proyecciones de State Plane, y Shapely no lo agobia con ellos.
Entonces lo combina con otros módulos de Python para leer o escribir archivos de forma y manipular proyecciones como osgeo.ogr, Fiona o PyShp .
Al buscar en Gis StackExchange, puede encontrar muchos ejemplos, pero le doy otro para ilustrar la combinación de shapely y Fiona y el uso de las funciones shapely intersection () y buffer () (Esto podría haberse hecho con PyShp).
Dados dos archivos de forma de polilínea:
Calcule la intersección (función intersección () de bien proporcionado)
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)Guarde el resultado como un nuevo archivo de forma
# 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}})Resultado:
Almacenar puntos individuales (función buffer () de bien proporcionado)
# 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}})Resultado
Almacenar la geometría de 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}})
Aquí está mi lista de software de geoprocesamiento Python.
Mi biblioteca 'ir a' de geoprocesamiento es la 'Biblioteca de sensores remotos y SIG' (RSGISLib). Es fácil de instalar y usar y la documentación es realmente buena. Tiene funcionalidad para el procesamiento de vectores y ráster: rara vez tengo que acercarme a una interfaz gráfica de usuario. Se puede encontrar aquí: http://rsgislib.org .
Un ejemplo en este caso es:
rsgislib.vectorutils.buffervector(inputvector, outputvector, bufferDist, force)Un comando para almacenar un vector en una memoria intermedia a una distancia especificada.
Dónde:
Ejemplo:
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)