4

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

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

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
4

2 回答 2

9

我喜欢这个问题。我怀疑我能给你最好的答案，而且绝对不能帮助 OGR，但是 FWIW 我会告诉你我现在在做什么。

我使用GeoPandas ，它是pandas的地理空间扩展。我推荐它——它是高级的，而且功能很多，免费为你提供Shapely和fiona的一切。twitter/@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)

开发人员告诉我，空间连接是“开发版本中的新功能”，所以如果你想在那里闲逛，我很想听听它是怎么回事！我的代码的主要问题是它很慢。

于 2014-12-22T19:11:21.123 回答
0 
import geopandas as gpd

# Read the data.
polygons = gpd.GeoDataFrame.from_file('polygons.shp')
points = gpd.GeoDataFrame.from_file('points.shp')

# Spatial Joins
pointsInPolygon = gpd.sjoin(points, polygons, how="inner", op='intersects')

# Add a field with 1 as a constant value
pointsInPolygon['const']=1

# Group according to the column by which you want to aggregate data
pointsInPolygon.groupby(['statename']).sum()
 
**The column ['const'] will give you the count number of points in your multipolygons.**

#If you want to see others columns as well, just type something like this : 
pointsInPolygon = pointsInPolygon.groupby('statename').agg({'columnA':'first', 'columnB':'first', 'const':'sum'}).reset_index()
    
    
  [1]: https://geopandas.org/docs/user_guide/mergingdata.html#spatial-joins
  [2]: https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.groupby.html
于 2021-03-07T16:36:05.393 回答