title: Python shapely: Aggregating points to shape files for a Choropleth map
body: <p>I'm trying a create a Choropleth in Python3 using shapely, fiona &amp; bokeh for display.</p>
  
  <p>I have a file with about 7000 lines that have the location of a town and a counter.</p>
  
  <p>Example:</p>
  
  <pre><code>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
  ...
  </code></pre>
  
  <p>I want to show these in a 12,5km grid, for which a shapefile is available on 
  <a href="https://opendata-esri-de.opendata.arcgis.com/datasets/3c1f46241cbb4b669e18b002e4893711_0" rel="nofollow noreferrer">https://opendata-esri-de.opendata.arcgis.com/datasets/3c1f46241cbb4b669e18b002e4893711_0</a></p>
  
  <p>The code I have works. </p>
  
  <p>It's very slow, because it's a brute force algorithm that checks each of the 7127 grid points against all of the 7000 points.</p>
  
  <pre><code>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]
  </code></pre>
  
  <p>Since this is extremly slow, I'm wondering if there is a faster way to get the <em>num_hex</em> value, especially, since the real world list of points may be a lot larger and a smaller grid with more shapes would deliver a better result.</p>