使用 Python 计算多多边形 shapefile 中的点数

2022-01-14 00:00:00 python shapefile shapely geolocation qgis

问题描述

我有一个由各个州组成的美国多边形 shapefile 作为它们的属性值.此外，我有存储我也感兴趣的点事件的纬度和经度值的数组.本质上，我想空间连接"点和多边形(或执行检查以查看每个多边形 [即状态]点在)，然后将每个状态的点数相加，找出哪个状态的事件"数量最多.

I have a polygon shapefile of the U.S. made up of individual states as their attribute values. In addition, I have arrays storing latitude and longitude values of point events that I am also interested in. Essentially, I would like to 'spatial join' the points and polygons (or perform a check to see which polygon [i.e., state] each point is in), then sum the number of points in each state to find out which state has the most number of 'events'.

我相信伪代码会是这样的:

I believe the pseudocode would be something like:

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: '

任何库或语法将不胜感激.

Any libraries or syntax would be greatly appreciated.

据我所知，OGR 库是我所需要的，但我在语法上有问题:

Based on what I can tell, the OGR library is what I need, but I am having trouble with the syntax:

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 like the question. I doubt I can give you the best answer, and definitely can't help with OGR, but FWIW I'll tell you what I'm doing right now.

我使用 GeoPandas，这是 熊猫.我推荐它——它是高级的并且功能很多，为您提供 Shapely 和 fiona 免费.twitter/@kajord 和其他人正在积极开发它.

I use GeoPandas, a geospatial extension of pandas. I recommend it — it's high-level and does a lot, giving you everything in Shapely and fiona for free. It is in active development by twitter/@kajord and others.

这是我的工作代码的一个版本.它假定您在 shapefile 中拥有所有内容，但很容易从列表中生成 geopandas.GeoDataFrame.

Here's a version of my working code. It assumes you have everything in shapefiles, but it's easy to generate a geopandas.GeoDataFrame from a list.

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)

开发人员告诉我，空间连接是开发版本中的新功能"，所以如果您想在 那里，我很想听听进展如何！我的代码的主要问题是它很慢.

The developer tells me that spatial joins are 'new in the dev version', so if you feel like poking around in there, I'd love to hear how that goes! The main problem with my code is that it's slow.

相关文章