À l'exception d'ArcPy, existe-t-il une bibliothèque python capable de faire du géotraitement, comme un tampon / intersection, avec des fichiers de formes?
16
À l'exception d'ArcPy, existe-t-il une bibliothèque python capable de faire du géotraitement, comme un tampon / intersection, avec des fichiers de formes?
Le livre de recettes Python GDAL / OGR contient un exemple de code pour mettre en mémoire tampon une géométrie .
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())et pour calculer l'intersection entre deux géométries
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()Les géométries peuvent être lues et écrites dans des fichiers de formes et une variété d'autres formats.
Pour simplifier, Shapely: manual permet tout le traitement géométrique de PostGIS en Python.
La première prémisse de Shapely est que les programmeurs Python devraient être capables d'effectuer des opérations de géométrie de type PostGIS en dehors d'un SGBDR ...
Le premier exemple 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), ....}L'exemple du polygone 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 deuxième prémisse est que la persistance, la sérialisation et la projection cartographique des entités sont des problèmes importants mais orthogonaux. Vous n'avez peut-être pas besoin d'une centaine de lecteurs et d'écrivains au format SIG ou de la multitude de projections State Plane, et Shapely ne vous en charge pas.
Vous le combinez donc avec d'autres modules Python pour lire ou écrire des fichiers de formes et manipuler des projections comme osgeo.ogr, Fiona ou PyShp .
En recherchant dans Gis StackExchange, vous pouvez trouver de nombreux exemples mais je vous en donne un autre pour illustrer la combinaison de shapely et Fiona et l'utilisation des fonctions shapely intersection () et buffer () (Cela aurait pu être fait avec PyShp).
Étant donné deux fichiers de formes polylignes:
Calculer l'intersection (fonction intersection () de galbée)
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)Enregistrez le résultat en tant que nouveau fichier de formes
# 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}})Résultat:
Tampon de points individuels (fonction buffer () de galbée)
# 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}})Résultat
Mettre en mémoire tampon la géométrie 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}})
Voici ma liste de logiciels de géotraitement Python.
Ma bibliothèque de géotraitement «go to» est la «bibliothèque de télédétection et SIG» (RSGISLib). C'est facile à installer et à utiliser et la documentation est vraiment bonne. Il a des fonctionnalités pour le traitement vectoriel et raster - je dois très rarement aller près d'un gui. Il peut être trouvé ici: http://rsgislib.org .
Un exemple dans ce cas est:
rsgislib.vectorutils.buffervector(inputvector, outputvector, bufferDist, force)Une commande pour tamponner un vecteur d'une distance spécifiée.
Où:
Exemple:
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)