1

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 ?

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.

  • 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 Jan 13 '17 at 23:13
  • can i use cities_with_country = gpd.sjoin(cities, countries, how="inner", op='intersects') where the cities and countries is the shapefiles ? – Chris Papas Jan 13 '17 at 23:15
  • 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 ? – Chris Papas Jan 15 '17 at 10:12
0

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')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.