我有一个美国的多边形形状文件,由各个州组成作为它们的属性值。此外,我有数组存储我也感兴趣的点事件的纬度和经度值。本质上,我想“空间连接”点和多边形(或执行检查以查看每个多边形[即状态]点在),然后将每个状态中的点数相加,以找出哪个状态具有最多的“事件”数量。
我相信伪代码会是这样的:
Read in US.shp
Read in lat/lon points of events
Loop through each state in the shapefile and find number of points in each state
print 'Here is a list of the number of points in each state: '
任何库或语法将不胜感激。
据我所知,OGR 库是我所需要的,但我在语法上遇到了问题:
dsPolygons = ogr.Open('US.shp')
polygonsLayer = dsPolygons.GetLayer()
#Iterating all the polygons
polygonFeature = polygonsLayer.GetNextFeature()
k=0
while polygonFeature:
k = k + 1
print "processing " + polygonFeature.GetField("STATE") + "-" + str(k) + " of " + str(polygonsLayer.GetFeatureCount())
geometry = polygonFeature.GetGeometryRef()
#Read in some points?
geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(-122.33,47.09)
point.AddPoint(-110.11,33.33)
#geomcol.AddGeometry(point)
print point.ExportToWkt()
print point
numCounts=0.0
while pointFeature:
if pointFeature.GetGeometryRef().Within(geometry):
numCounts = numCounts + 1
pointFeature = pointsLayer.GetNextFeature()
polygonFeature = polygonsLayer.GetNextFeature()
#Loop through to see how many events in each state
我喜欢这个问题。我怀疑我能给你最好的答案,而且绝对不能帮助 OGR,但 FWIW 我会告诉你我现在在做什么。
I use 地理熊猫,地理空间延伸pandas。我推荐它——它的水平很高,功能很多,可以为您提供一切Shapely and fiona免费。它正在积极开发中推特/@kajord和别的。
这是我的工作代码的一个版本。它假设您拥有 shapefile 中的所有内容,但很容易生成geopandas.GeoDataFrame
从列表中。
import geopandas as gpd
# Read the data.
polygons = gpd.GeoDataFrame.from_file('polygons.shp')
points = gpd.GeoDataFrame.from_file('points.shp')
# Make a copy because I'm going to drop points as I
# assign them to polys, to speed up subsequent search.
pts = points.copy()
# We're going to keep a list of how many points we find.
pts_in_polys = []
# Loop over polygons with index i.
for i, poly in polygons.iterrows():
# Keep a list of points in this poly
pts_in_this_poly = []
# Now loop over all points with index j.
for j, pt in pts.iterrows():
if poly.geometry.contains(pt.geometry):
# Then it's a hit! Add it to the list,
# and drop it so we have less hunting.
pts_in_this_poly.append(pt.geometry)
pts = pts.drop([j])
# We could do all sorts, like grab a property of the
# points, but let's just append the number of them.
pts_in_polys.append(len(pts_in_this_poly))
# Add the number of points for each poly to the dataframe.
polygons['number of points'] = gpd.GeoSeries(pts_in_polys)
开发人员告诉我,空间连接是“开发版本中的新功能”,所以如果您想四处看看in there,我很想听听事情进展如何!我的代码的主要问题是它很慢。
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)