Python을 사용하여 다중 다각형 모양 파일의 점 수를 계산합니다.

Question

미국의 다각형 모양 파일은 특성 상태 값으로 개별 상태로 구성되어 있습니다. 또한, 관심이있는 지점 이벤트의 위도와 경도 값을 저장하는 배열이 있습니다. 본질적으로, 점과 다각형을 '공간 결합'하고 싶습니다 (또는 각 폴리곤 [즉, 상태]를 확인하는 검사를 수행합니다 포인트가있는 경우), 각 상태의 포인트 수를 합하여 가장 많은 수의 '이벤트'가있는 상태를 찾습니다.

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 

Answer 1

내가 질문을 좋아한다. 나는 최선의 대답을 줄 수 있을지 의심 스럽지만 확실히 OGR을 도울 수는 없지만 FWIW는 내가 지금하고있는 일을 말해 줄 것이다.

나는 pandas의 지형 공간 확장 인 GeoPandas을 사용합니다. 나는 그것을 추천한다 - 그것은 높은 수준이고, 많은 것을한다, 당신에게 Shapely와 fiona의 모든 것을 무료로 준다. 그것은 twitter/@kajord 및 다른 사람에 의해 활동적인 발달에있다.

여기 내 작업 코드 버전입니다. 그것은 당신이 shapefiles에있는 모든 것을 가지고 있다고 가정하지만 목록에서 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을 파고 같은 느낌, 그래서 만약 'dev에 버전의 새로운'되는 공간 조인 알려줍니다, 나는 어떻게되는지 듣고 싶어요! 내 코드의 주요 문제점은 느리다는 것이다.