Spatial join using python shapely and fiona

Question

I want to create a python script when count how many points(from point shapefile) is in the polygons(polygon shapefile) and I want that count to create a new polygon shapefiles with the polygons where have upper to one points within.

using fiona for inputs shapefiles

filepath1 = "polygon.shp"
file1features = []
outschema = None
with fiona.collection(filepath1, "r") as input1:
 outschema = input1.schema.copy()
 for p1 in input1:
 file1features.append(p1)
filepath2 = "point.shp"
file2features = []
outschema = None
with fiona.collection(filepath2, "r") as input1:
 outschema = input2.schema.copy()
 for p1 in input2:
 file2features.append(p1)

I have find this code from github for spatial join but how to connect that code woth my code for inputs shapefiles and how to crate a new polygon shapefile with upper count> 1 ?

import numpy as np
import pandas as pd
from shapely import prepared
def sjoin(left_df, right_df, how='inner', op='intersects',
 lsuffix='left', rsuffix='right'):
 """Spatial join of two GeoDataFrames.
 Parameters
 ----------
 left_df, right_df : GeoDataFrames
 how : string, default 'inner'
 The type of join:
 * 'left': use keys from left_df; retain only left_df geometry column
 * 'right': use keys from right_df; retain only right_df geometry column
 * 'inner': use intersection of keys from both dfs; retain only
 left_df geometry column
 op : string, default 'intersection'
 Binary predicate, one of {'intersects', 'contains', 'within'}.
 See http://toblerity.org/shapely/manual.html#binary-predicates.
 lsuffix : string, default 'left'
 Suffix to apply to overlapping column names (left GeoDataFrame).
 rsuffix : string, default 'right'
 Suffix to apply to overlapping column names (right GeoDataFrame).
 """
 import rtree
 allowed_hows = ['left', 'right', 'inner']
 if how not in allowed_hows:
 raise ValueError("`how` was \"%s\" but is expected to be in %s" % \
 (how, allowed_hows))
 allowed_ops = ['contains', 'within', 'intersects']
 if op not in allowed_ops:
 raise ValueError("`op` was \"%s\" but is expected to be in %s" % \
 (op, allowed_ops))
 if op == "within":
 # within implemented as the inverse of contains; swap names
 left_df, right_df = right_df, left_df
 if left_df.crs != right_df.crs:
 print('Warning: CRS does not match!')
 tree_idx = rtree.index.Index()
 right_df_bounds = right_df.geometry.apply(lambda x: x.bounds)
 for i in right_df_bounds.index:
 tree_idx.insert(i, right_df_bounds[i])
 idxmatch = (left_df.geometry.apply(lambda x: x.bounds)
 .apply(lambda x: list(tree_idx.intersection(x))))
 idxmatch = idxmatch[idxmatch.apply(len) > 0]
 if idxmatch.shape[0] > 0:
 # if output from join has overlapping geometries
 r_idx = np.concatenate(idxmatch.values)
 l_idx = np.concatenate([[i] * len(v) for i, v in idxmatch.iteritems()])
 # Vectorize predicate operations
 def find_intersects(a1, a2):
 return a1.intersects(a2)
 def find_contains(a1, a2):
 return a1.contains(a2)
 predicate_d = {'intersects': find_intersects,
 'contains': find_contains,
 'within': find_contains}
 check_predicates = np.vectorize(predicate_d[op])
 result = (
 pd.DataFrame(
 np.column_stack(
 [l_idx,
 r_idx,
 check_predicates(
 left_df.geometry
 .apply(lambda x: prepared.prep(x))[l_idx],
 right_df[right_df.geometry.name][r_idx])
 ]))
 )
 result.columns = ['index_%s' % lsuffix, 'index_%s' % rsuffix, 'match_bool']
 result = (
 pd.DataFrame(result[result['match_bool']==1])
 .drop('match_bool', axis=1)
 )
 else:
 # when output from the join has no overlapping geometries
 result = pd.DataFrame(columns=['index_%s' % lsuffix, 'index_%s' % rsuffix])
 if op == "within":
 # within implemented as the inverse of contains; swap names
 left_df, right_df = right_df, left_df
 result = result.rename(columns={
 'index_%s' % (lsuffix): 'index_%s' % (rsuffix),
 'index_%s' % (rsuffix): 'index_%s' % (lsuffix)})
 if how == 'inner':
 result = result.set_index('index_%s' % lsuffix)
 return (
 left_df
 .merge(result, left_index=True, right_index=True)
 .merge(right_df.drop(right_df.geometry.name, axis=1),
 left_on='index_%s' % rsuffix, right_index=True,
 suffixes=('_%s' % lsuffix, '_%s' % rsuffix))
 )
 elif how == 'left':
 result = result.set_index('index_%s' % lsuffix)
 return (
 left_df
 .merge(result, left_index=True, right_index=True, how='left')
 .merge(right_df.drop(right_df.geometry.name, axis=1),
 how='left', left_on='index_%s' % rsuffix, right_index=True,
 suffixes=('_%s' % lsuffix, '_%s' % rsuffix))
 )
 elif how == 'right':
 return (
 left_df
 .drop(left_df.geometry.name, axis=1)
 .merge(result.merge(right_df,
 left_on='index_%s' % rsuffix, right_index=True,
 how='right'), left_index=True,
 right_on='index_%s' % lsuffix, how='right')
 .set_index('index_%s' % rsuffix)
)

I thing so like this for new shapefile ?

if intersectfeatures:
 outfile = "utputfile.shp"
 with fiona.collection(outfile, "w", "ESRI Shapefile", outschema) as output:
 for outfeature in intersectfeatures:
 output.write(outfeature)

But I can't to connection, any idea ?

Answer 1

For the first, you could to create a spatial join where you merge points to polygons, so you'd get a frequency of points per polygon as an attribute. Spatial joins are covered at geopandas here and are reviewed in a simple notebook here. The code snippet you include from their site is more complicated, and what's "under the hood" of their library. You could add geopandas as a dependency, and that will even reduce your fiona code.

For the second goal you state, do you really need a new shapefile with counts greater than 0, or is it a visualization task (where you select to visualize polygons with 1 or more points)? If you need a new shapefile, you'd select your polygons with frequency> 0, save as a GeoDataFrame, and then export as a shapefile.

Note: You'll need to also install rtree to get the spatial join to work properly in geopandas.

• Chris Papas: yes i need new shapefile but for the 'select your polygons with frequency > 0, save as a GeoDataFrame' must do it after from spatial join right?i dont care the progress i need fast calculator
• Chris Papas: can i use cities_with_country = gpd.sjoin(cities, countries, how="inner", op='intersects') where the cities and countries is the shapefiles ?
• Chris Papas: i thing so i find it with geopandas you are right is a very simple on my task i thing so i need merge but how to count the points in the polygon ?

Answer 2

Try this code, it should work

import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
from fiona.crs import from_epsg
import fiona
df = gpd.sjoin(points, polys[['column1', 'column2', 'geometry']], how='left', op='intersects') # Here I am using intersects with left join

If you are using op='within' then not use left like

df = gpd.sjoin(points, polys[['column1', 'column2', 'geometry']], op='within')