Python shapely: агрегирование точек в файлы форм для карты Choropleth - PullRequest
Новая
       48

Python shapely: агрегирование точек в файлы форм для карты Choropleth

0 голосов
/ 17 апреля 2019

Я пытаюсь создать Choropleth в Python3, используя shapely, fiona & bokeh для отображения.

У меня есть файл с примерно 7000 строками, в котором указано местоположение города и счетчик.

Пример: 

54.7604;9.55827;208
54.4004;9.95918;207
53.8434;9.95271;203
53.5979;10.0013;201
53.728;10.2526;197
53.646;10.0403;196
54.3977;10.1054;193
52.4385;9.39217;193
53.815;10.3476;192
...

Я хочу показать их в сетке 12,5 км, для которой доступен шейп-файл на https://opendata -esri-de.opendata.arcgis.com / наборы данных / 3c1f46241cbb4b669e18b002e4893711_0

Код у меня работает.

Это очень медленно, потому что это алгоритм грубой силы, который проверяет каждую из 7127 точек сетки на все 7000 точек. 

import pandas as pd
import fiona
from shapely.geometry import Polygon, Point, MultiPoint, MultiPolygon
from shapely.prepared import prep
sf = r'c:\Temp\geo_de\Hexagone_125_km\Hexagone_125_km.shp'
shp = fiona.open(sf)

district_xy = [ [ xy for xy in feat["geometry"]["coordinates"][0]] for feat in shp] 
district_poly = [ Polygon(xy) for xy in district_xy] # coords to Polygon

df_p = pd.read_csv('points_file.csv', sep=';', header=None)
df_p.columns = ('lat', 'lon', 'count')

map_points = [Point(x,y) for x,y in zip(df_p.lon, df_p.lat)] # Convert Points to Shapely Points

all_points = MultiPoint(map_points) # all points

def calc_points_per_poly(poly, points, values): # Returns total for poly
    poly = prep(poly)
    return sum([v for p, v in zip(points, values) if poly.contains(p)])

# this is the slow part
# for each shape this sums um the points

sum_hex = [calc_points_per_poly(x, all_points, df_p['count']) for x in district_poly]

Так как это очень медленно, мне интересно, есть ли более быстрый способ получить значение num_hex , особенно, поскольку список точек реального мира может быть намного больше и меньше сетка с больше форм дало бы лучший результат.

1 Ответ

1 голос
/ 18 апреля 2019

Я бы порекомендовал использовать «геопанда» и его встроенный пространственный индекс rtree.Это позволяет вам выполнять проверку только в том случае, если существует вероятность того, что точка лежит в пределах многоугольника. 

import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon, Point

sf = 'Hexagone_125_km.shp'
shp = gpd.read_file(sf)

df_p = pd.read_csv('points_file.csv', sep=';', header=None)
df_p.columns = ('lat', 'lon', 'count')

gdf_p = gpd.GeoDataFrame(df_p, geometry=[Point(x,y) for x,y in zip(df_p.lon, df_p.lat)])

sum_hex = []
spatial_index = gdf_p.sindex

for index, row in shp.iterrows():
    polygon = row.geometry
    possible_matches_index = list(spatial_index.intersection(polygon.bounds))
    possible_matches = gdf_p.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.within(polygon)]
    sum_hex.append(sum(precise_matches['count']))

shp['sum'] = sum_hex

Это решение должно быть быстрее вашего.Затем вы можете построить свой GeoDataFrame через Bokeh.Если вы хотите узнать больше о пространственной индексации, я рекомендую эту статью Джеффа Боинга: https://geoffboeing.com/2016/10/r-tree-spatial-index-python/

...