矢量操作可以使用 ogr这类基础方法进行,但同时还有一些高级库支持更便捷的矢量操作,例如Fiona(shapely)、pyshp(shapefile)、geopandas等。
教程资料导航
ogr教程资料:
1.GDAL/OGR例子合集
2.GDAL官网API
3.网上的学习笔记,教程来自于犹他州立大学
4.中国OSGeo教程
shapely教程资料:
pyshp教程资料:
GeoPandas教程资料:
Fiona教程资料:
以下为判断点在面内的示例:
Solution 1 with Fiona
import fiona
from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
multipol = fiona.open("multipol.shp")
multi= multipol.next() # only one feature in the shapefile
print multi
{'geometry': {'type': 'MultiPolygon', 'coordinates': [[[(-0.5275288092189501, 0.5569782330345711), (-0.11779769526248396, 0.29065300896286816), (-0.25608194622279135, 0.01920614596670933), (-0.709346991037132, -0.08834827144686286), (-0.8629961587708066, 0.18309859154929575), (-0.734955185659411, 0.39820742637644047), (-0.5275288092189501, 0.5569782330345711)]], [[(0.19974391805377723, 0.060179257362355965), (0.5480153649167734, 0.1293213828425096), (0.729833546734955, 0.03969270166453265), (0.8143405889884763, -0.13956466069142115), (0.701664532650448, -0.38540332906530095), (0.4763124199743918, -0.5006402048655569), (0.26888604353393086, -0.4238156209987196), (0.18950064020486557, -0.2291933418693981), (0.19974391805377723, 0.060179257362355965)]], [[(-0.3764404609475033, -0.295774647887324), (-0.11523687580025621, -0.3597951344430217), (-0.033290653008962945, -0.5800256081946222), (-0.11523687580025621, -0.7413572343149808), (-0.3072983354673495, -0.8591549295774648), (-0.58898847631242, -0.6927016645326505), (-0.6555697823303457, -0.4750320102432779), (-0.3764404609475033, -0.295774647887324)]]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}
Fiona interprets the feature as a MultiPolygon and you can apply the solution presented in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (1)
points= ([pt for pt in fiona.open("points.shp")])
for i, pt in enumerate(points):
point = shape(pt['geometry'])
if point.within(shape(multi['geometry'])):
print i, shape(points[i]['geometry'])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)
Solution 2 with pyshp (shapefile) and the geo_interface (GeoJSON like) protocol
import shapefile
pts = shapefile.Reader("points.shp")
polys = shapefile.Reader("multipol.shp")
points = [pt.shape.__geo_interface__ for pt in pts.shapeRecords()]
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__) # 1 polygon
print multi
MULTIPOLYGON (((-0.5275288092189501 0.5569782330345711, -0.117797695262484 0.2906530089628682, -0.2560819462227913 0.01920614596670933, -0.7093469910371319 -0.08834827144686286, -0.8629961587708066 0.1830985915492958, -0.734955185659411 0.3982074263764405, -0.5275288092189501 0.5569782330345711)), ((0.1997439180537772 0.06017925736235596, 0.5480153649167734 0.1293213828425096, 0.729833546734955 0.03969270166453265, 0.8143405889884763 -0.1395646606914211, 0.701664532650448 -0.3854033290653009, 0.4763124199743918 -0.5006402048655569, 0.2688860435339309 -0.4238156209987196, 0.1895006402048656 -0.2291933418693981, 0.1997439180537772 0.06017925736235596)), ((-0.3764404609475033 -0.295774647887324, -0.1152368758002562 -0.3597951344430217, -0.03329065300896294 -0.5800256081946222, -0.1152368758002562 -0.7413572343149808, -0.3072983354673495 -0.8591549295774648, -0.58898847631242 -0.6927016645326505, -0.6555697823303457 -0.4750320102432779, -0.3764404609475033 -0.295774647887324)))
for i, pt in enumerate(points):
point = shape(pt)
if point.within(multi):
print i, shape(points[i])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)
Solution 3 with ogr and the geo_interface protocol (Python Geo_interface applications)
from osgeo import ogr
import json
def records(file):
# generator
reader = ogr.Open(file)
layer = reader.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
yield json.loads(feature.ExportToJson())
points = [pt for pt in records("point_multi_contains.shp")]
multipol = records("multipol.shp")
multi = multipol.next() # 1 feature
for i, pt in enumerate(points):
point = shape(pt['geometry'])
if point.within(shape(multi['geometry'])):
print i, shape(points[i]['geometry'])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.499359795134443 -0.060179257362356)
5 POINT (-0.376440460947503 -0.475032010243278)
6 POINT (-0.309859154929577 -0.631241997439181)
Solution 4 with GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)
import geopandas
point = geopandas.GeoDataFrame.from_file('points.shp')
poly = geopandas.GeoDataFrame.from_file('multipol.shp')
from geopandas.tools import sjoin
pointInPolys = sjoin(point, poly, how='left')
grouped = pointInPolys.groupby('index_right')
list(grouped)
[(0.0, geometry id_left index_right id_right
1 POINT (-0.58898847631242 0.17797695262484) None 0.0 1.0
3 POINT (0.4993597951344431 -0.06017925736235585) None 0.0 1.0
5 POINT (-0.3764404609475033 -0.4750320102432779) None 0.0 1.0
6 POINT (-0.3098591549295775 -0.6312419974391805) None 0.0 1.0 ]
print grouped.groups
{0.0: [1, 3, 5, 6]}
The points 1,3,5,6 falls within the boundaries of the MultiPolygon
参考:https://gis.stackovernet.com/cn/q/56553