From 33e9aa5dc7ddff50c87506ad5e7fe00ca57f7626 Mon Sep 17 00:00:00 2001 From: mdbartos Date: Sat, 4 Apr 2015 03:38:51 -0700 Subject: [PATCH 1/8] Optimize sjoin --- geopandas/tools/sjoin.py | 92 +++++++++++++++++++++++++--------------- 1 file changed, 58 insertions(+), 34 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index 896d823..634e465 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -1,8 +1,10 @@ +import geopandas as gpd +import numpy as np import pandas as pd -from geopandas import GeoDataFrame -from .overlay import _uniquify +import rtree +from shapely import prepared -def sjoin(left_df, right_df, how="left", op="intersects", use_sindex=True, **kwargs): +def sjoin(left_df, right_df, how='left', op='intersects', crs_convert=True, lsuffix='left', rsuffix='right', **kwargs): """Spatial join of two GeoDataFrames. left_df, right_df are GeoDataFrames @@ -15,48 +17,70 @@ def sjoin(left_df, right_df, how="left", op="intersects", use_sindex=True, **kwa use_sindex : Use the spatial index to speed up operation? Default is True kwargs: passed to op method """ + + # CHECK VALIDITY OF JOIN TYPE 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)) + + # CHECK VALIDITY OF PREDICATE OPERATION + allowed_ops = ['contains', 'within', 'intersects'] - if how == "right": - # right outer join just implemented as the inverse of left; swap names + if op not in allowed_ops: + raise ValueError("`op` was \"%s\" but is expected to be in %s" % \ + (op, allowed_ops)) + + # IF WITHIN, SWAP NAMES + if op == "within": + # within implemented as the inverse of contains; swap names left_df, right_df = right_df, left_df - collection = [] - for i, feat in left_df.iterrows(): - geom = feat.geometry + # CONVERT CRS IF NOT EQUAL + if left_df.crs != right_df.crs: + print 'Warning: CRS does not match!' + if crs_convert == True: + print 'Converting CRS...' + if left_df.values.nbytes >= right_df.values.nbytes: + right_df = right_df.to_crs(left_df.crs) + elif left_df.values.nbytes < right_df.values.nbytes: + left_df = left_df.to_crs(right_df.crs) - if use_sindex and right_df.sindex: - candidates = [x.object for x in - right_df.sindex.intersection(geom.bounds, objects=True)] - else: - candidates = [i for i, x in right_df.iterrows()] + # CONSTRUCT SPATIAL INDEX FOR RIGHT DATAFRAME + 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]) - feature_hits = 0 - for cand_id in candidates: - candidate = right_df.ix[cand_id] - if getattr(geom, op)(candidate.geometry, **kwargs): - newseries = candidate.drop(right_df._geometry_column_name) - newfeat = pd.concat([feat, newseries]) - newfeat.index = _uniquify(newfeat.index) - collection.append(newfeat) - feature_hits += 1 + # FIND INTERSECTION OF SPATIAL INDEX + idxmatch = left_df['geometry'].apply(lambda x: x.bounds).apply(lambda x: list(tree_idx.intersection(x))) + idxmatch = idxmatch[idxmatch.str.len() > 0] - # TODO Should we perform aggregation if feature_hit > 1? - # Advantage: single step and possible performance improvement - # Disadvantage: Pandas already has groupby so user can do this later + r_idx = np.concatenate(idxmatch.values) + l_idx = np.concatenate((idxmatch.str.len()*pd.Series([[i] for i in idxmatch.index], index=idxmatch.index)).values) - # If left does not spatially join with any right features, - # Fill in the right columns with NA - if how != 'inner' and feature_hits == 0: - empty = pd.Series(dict.fromkeys(right_df.columns, None)) - empty.drop(right_df._geometry_column_name, inplace=True) + # VECTORIZE PREDICATE OPERATIONS + def find_intersects(a1, a2): + return a1.intersects(a2) - newfeat = pd.concat([feat, empty]) - newfeat.index = _uniquify(newfeat.index) - collection.append(newfeat) + def find_contains(a1, a2): + return a1.contains(a2) - return GeoDataFrame(collection, index=range(len(collection))) + predicate_d = {'intersects': find_intersects, 'contains': find_contains, 'within': find_contains} + + check_predicates = np.vectorize(predicate_d[op]) + + # CHECK PREDICATES + result = pd.DataFrame(np.column_stack([l_idx, r_idx, check_predicates(left_df['geometry'].apply(lambda x: prepared.prep(x)).values[l_idx], right_df['geometry'].values[r_idx])])) + result.columns = ['index_%s' % lsuffix, 'index_%s' % rsuffix, 'match_bool'] + result = pd.DataFrame(result[result['match_bool']==1].set_index('index_%s' % lsuffix)['index_%s' % rsuffix]) + + # IF 'WITHIN', SWAP NAMES AGAIN + if op == "within": + # within implemented as the inverse of contains; swap names + left_df, right_df = right_df, left_df + result = result.reset_index().rename(columns={'index_%s' % (lsuffix): 'index_%s' % (rsuffix), 'index_%s' % (rsuffix): 'index_%s' % (lsuffix)}).set_index('index_left').sort_index() + + # APPLY JOIN + return left_df.merge(result, left_index=True, right_index=True).merge(right_df, left_on='index_%s' % rsuffix, right_index=True, how=how, suffixes=('_%s' % lsuffix, '_%s' % rsuffix)) From 17c67d13a263fdca8bf0f45a609e65b9142b5da6 Mon Sep 17 00:00:00 2001 From: mdbartos Date: Sat, 4 Apr 2015 05:15:28 -0700 Subject: [PATCH 2/8] Marked test file for sources of error --- geopandas/tools/sjoin.py | 2 ++ tests/test_sjoin.py | 49 ++++++++++++++++++++-------------------- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index 634e465..e03f15b 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -46,6 +46,8 @@ def sjoin(left_df, right_df, how='left', op='intersects', crs_convert=True, lsuf right_df = right_df.to_crs(left_df.crs) elif left_df.values.nbytes < right_df.values.nbytes: left_df = left_df.to_crs(right_df.crs) + else: + pass # CONSTRUCT SPATIAL INDEX FOR RIGHT DATAFRAME tree_idx = rtree.index.Index() diff --git a/tests/test_sjoin.py b/tests/test_sjoin.py index 219d7ab..69db7b0 100644 --- a/tests/test_sjoin.py +++ b/tests/test_sjoin.py @@ -25,51 +25,52 @@ class TestSpatialJoin(unittest.TestCase): shutil.rmtree(self.tempdir) def test_sjoin_left(self): - df = sjoin(self.pointdf, self.polydf) - self.assertEquals(df.shape, (21,7)) - for i, row in df.iterrows(): - self.assertEquals(row.geometry.type, 'Point') + df = sjoin(self.pointdf, self.polydf, crs_convert=False) + self.assertEquals(df.shape, (11,9)) +# for i, row in df.iterrows(): +# self.assertEquals(row.geometry.type, 'Point') self.assertTrue('pointattr1' in df.columns) self.assertTrue('BoroCode' in df.columns) def test_sjoin_right(self): # the inverse of left - df = sjoin(self.pointdf, self.polydf, how="right") - df2 = sjoin(self.polydf, self.pointdf, how="left") - self.assertEquals(df.shape, (12, 7)) - self.assertEquals(df.shape, df2.shape) - for i, row in df.iterrows(): - self.assertEquals(row.geometry.type, 'MultiPolygon') - for i, row in df2.iterrows(): - self.assertEquals(row.geometry.type, 'MultiPolygon') + df = sjoin(self.pointdf, self.polydf, how="right", crs_convert=False) + df2 = sjoin(self.polydf, self.pointdf, how="left", crs_convert=False) + self.assertEquals(df.shape, (12, 9)) +# self.assertEquals(df.shape, df2.shape) +# for i, row in df.iterrows(): +# self.assertEquals(row.geometry.type, 'MultiPolygon') +# for i, row in df2.iterrows(): +# self.assertEquals(row.geometry.type, 'MultiPolygon') def test_sjoin_inner(self): - df = sjoin(self.pointdf, self.polydf, how="inner") - self.assertEquals(df.shape, (11, 7)) + df = sjoin(self.pointdf, self.polydf, how="inner", crs_convert=False) + self.assertEquals(df.shape, (11, 9)) def test_sjoin_op(self): # points within polygons - df = sjoin(self.pointdf, self.polydf, how="left", op="within") - self.assertEquals(df.shape, (21,7)) + df = sjoin(self.pointdf, self.polydf, how="left", op="within", crs_convert=False) + self.assertEquals(df.shape, (11,9)) self.assertAlmostEquals(df.ix[1]['Shape_Leng'], 330454.175933) # points contain polygons? never happens so we should have nulls - df = sjoin(self.pointdf, self.polydf, how="left", op="contains") - self.assertEquals(df.shape, (21, 7)) - self.assertEquals(df.ix[1]['Shape_Area'], None) +# df = sjoin(self.pointdf, self.polydf, how="left", op="contains", crs_convert=False) +# self.assertEquals(df.shape, (11, 9)) +# self.assertEquals(df.ix[1]['Shape_Area'], None) - def test_sjoin_bad_op(self): + def test_sjoin_bad_op(self, crs_convert=False): # AttributeError: 'Point' object has no attribute 'spandex' - self.assertRaises(AttributeError, sjoin, + self.assertRaises(ValueError, sjoin, self.pointdf, self.polydf, how="left", op="spandex") - def test_sjoin_duplicate_column_name(self): + @unittest.skip("Not implemented") + def test_sjoin_duplicate_column_name(self, crs_convert=False): pointdf2 = self.pointdf.rename(columns={'pointattr1': 'Shape_Area'}) - df = sjoin(pointdf2, self.polydf, how="left") + df = sjoin(pointdf2, self.polydf, how="left", crs_convert=False) self.assertTrue('Shape_Area' in df.columns) self.assertTrue('Shape_Area_2' in df.columns) @unittest.skip("Not implemented") def test_sjoin_outer(self): df = sjoin(self.pointdf, self.polydf, how="outer") - self.assertEquals(df.shape, (21,7)) + self.assertEquals(df.shape, (21,9)) From dbe054ab0744d350d90fd77b24d4b545476bc5f7 Mon Sep 17 00:00:00 2001 From: mdbartos Date: Sat, 4 Apr 2015 17:38:26 -0700 Subject: [PATCH 3/8] Maintained compatibility with old tests. Edited sjoin to maintain python 3 compatibility. --- geopandas/tools/sjoin.py | 21 ++++++++++----- tests/test_sjoin.py | 57 ++++++++++++++++++++-------------------- 2 files changed, 43 insertions(+), 35 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index e03f15b..918bc34 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -4,7 +4,7 @@ import pandas as pd import rtree from shapely import prepared -def sjoin(left_df, right_df, how='left', op='intersects', crs_convert=True, lsuffix='left', rsuffix='right', **kwargs): +def sjoin(left_df, right_df, how='left', op='intersects', convert_crs=True, lsuffix='left', rsuffix='right', **kwargs): """Spatial join of two GeoDataFrames. left_df, right_df are GeoDataFrames @@ -39,9 +39,9 @@ def sjoin(left_df, right_df, how='left', op='intersects', crs_convert=True, lsuf # CONVERT CRS IF NOT EQUAL if left_df.crs != right_df.crs: - print 'Warning: CRS does not match!' - if crs_convert == True: - print 'Converting CRS...' + print('Warning: CRS does not match!') + if convert_crs == True: + print('Converting CRS...') if left_df.values.nbytes >= right_df.values.nbytes: right_df = right_df.to_crs(left_df.crs) elif left_df.values.nbytes < right_df.values.nbytes: @@ -76,13 +76,20 @@ def sjoin(left_df, right_df, how='left', op='intersects', crs_convert=True, lsuf # CHECK PREDICATES result = pd.DataFrame(np.column_stack([l_idx, r_idx, check_predicates(left_df['geometry'].apply(lambda x: prepared.prep(x)).values[l_idx], right_df['geometry'].values[r_idx])])) result.columns = ['index_%s' % lsuffix, 'index_%s' % rsuffix, 'match_bool'] - result = pd.DataFrame(result[result['match_bool']==1].set_index('index_%s' % lsuffix)['index_%s' % rsuffix]) + result = pd.DataFrame(result[result['match_bool']==1]).drop('match_bool', axis=1) # IF 'WITHIN', SWAP NAMES AGAIN if op == "within": # within implemented as the inverse of contains; swap names left_df, right_df = right_df, left_df - result = result.reset_index().rename(columns={'index_%s' % (lsuffix): 'index_%s' % (rsuffix), 'index_%s' % (rsuffix): 'index_%s' % (lsuffix)}).set_index('index_left').sort_index() + result = result.rename(columns={'index_%s' % (lsuffix): 'index_%s' % (rsuffix), 'index_%s' % (rsuffix): 'index_%s' % (lsuffix)}) # APPLY JOIN - return left_df.merge(result, left_index=True, right_index=True).merge(right_df, left_on='index_%s' % rsuffix, right_index=True, how=how, suffixes=('_%s' % lsuffix, '_%s' % rsuffix)) + 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('geometry', 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('geometry', axis=1), how='left', left_on='index_%s' % rsuffix, right_index=True, suffixes=('_%s' % lsuffix, '_%s' % rsuffix)) + elif how == 'right': + return left_df.drop('geometry', 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) diff --git a/tests/test_sjoin.py b/tests/test_sjoin.py index 69db7b0..34c0d4a 100644 --- a/tests/test_sjoin.py +++ b/tests/test_sjoin.py @@ -1,6 +1,8 @@ + from __future__ import absolute_import import tempfile import shutil +import numpy as np from shapely.geometry import Point from geopandas import GeoDataFrame, read_file from geopandas.tools import sjoin @@ -25,52 +27,51 @@ class TestSpatialJoin(unittest.TestCase): shutil.rmtree(self.tempdir) def test_sjoin_left(self): - df = sjoin(self.pointdf, self.polydf, crs_convert=False) - self.assertEquals(df.shape, (11,9)) -# for i, row in df.iterrows(): -# self.assertEquals(row.geometry.type, 'Point') + df = sjoin(self.pointdf, self.polydf, convert_crs=False) + self.assertEquals(df.shape, (21,8)) + for i, row in df.iterrows(): + self.assertEquals(row.geometry.type, 'Point') self.assertTrue('pointattr1' in df.columns) self.assertTrue('BoroCode' in df.columns) def test_sjoin_right(self): # the inverse of left - df = sjoin(self.pointdf, self.polydf, how="right", crs_convert=False) - df2 = sjoin(self.polydf, self.pointdf, how="left", crs_convert=False) - self.assertEquals(df.shape, (12, 9)) -# self.assertEquals(df.shape, df2.shape) -# for i, row in df.iterrows(): -# self.assertEquals(row.geometry.type, 'MultiPolygon') -# for i, row in df2.iterrows(): -# self.assertEquals(row.geometry.type, 'MultiPolygon') + df = sjoin(self.pointdf, self.polydf, how="right", convert_crs=False) + df2 = sjoin(self.polydf, self.pointdf, how="left", convert_crs=False) + self.assertEquals(df.shape, (12, 8)) + self.assertEquals(df.shape, df2.shape) + for i, row in df.iterrows(): + self.assertEquals(row.geometry.type, 'MultiPolygon') + for i, row in df2.iterrows(): + self.assertEquals(row.geometry.type, 'MultiPolygon') def test_sjoin_inner(self): - df = sjoin(self.pointdf, self.polydf, how="inner", crs_convert=False) - self.assertEquals(df.shape, (11, 9)) + df = sjoin(self.pointdf, self.polydf, how="inner", convert_crs=False) + self.assertEquals(df.shape, (11, 8)) def test_sjoin_op(self): # points within polygons - df = sjoin(self.pointdf, self.polydf, how="left", op="within", crs_convert=False) - self.assertEquals(df.shape, (11,9)) + df = sjoin(self.pointdf, self.polydf, how="left", op="within", convert_crs=False) + self.assertEquals(df.shape, (21,8)) self.assertAlmostEquals(df.ix[1]['Shape_Leng'], 330454.175933) # points contain polygons? never happens so we should have nulls -# df = sjoin(self.pointdf, self.polydf, how="left", op="contains", crs_convert=False) -# self.assertEquals(df.shape, (11, 9)) -# self.assertEquals(df.ix[1]['Shape_Area'], None) + df = sjoin(self.pointdf, self.polydf, how="left", op="contains", convert_crs=False) + self.assertEquals(df.shape, (21, 8)) + self.assertTrue(np.isnan(df.ix[1]['Shape_Area'])) - def test_sjoin_bad_op(self, crs_convert=False): + def test_sjoin_bad_op(self): # AttributeError: 'Point' object has no attribute 'spandex' self.assertRaises(ValueError, sjoin, - self.pointdf, self.polydf, how="left", op="spandex") + self.pointdf, self.polydf, how="left", op="spandex", convert_crs=False) - @unittest.skip("Not implemented") - def test_sjoin_duplicate_column_name(self, crs_convert=False): + def test_sjoin_duplicate_column_name(self): pointdf2 = self.pointdf.rename(columns={'pointattr1': 'Shape_Area'}) - df = sjoin(pointdf2, self.polydf, how="left", crs_convert=False) - self.assertTrue('Shape_Area' in df.columns) - self.assertTrue('Shape_Area_2' in df.columns) + df = sjoin(pointdf2, self.polydf, how="left", convert_crs=False) + self.assertTrue('Shape_Area_left' in df.columns) + self.assertTrue('Shape_Area_right' in df.columns) @unittest.skip("Not implemented") def test_sjoin_outer(self): - df = sjoin(self.pointdf, self.polydf, how="outer") - self.assertEquals(df.shape, (21,9)) + df = sjoin(self.pointdf, self.polydf, how="outer", convert_crs=False) + self.assertEquals(df.shape, (21,8)) From 1cdad116792198c7804f76949139c9029f5df5b1 Mon Sep 17 00:00:00 2001 From: mdbartos Date: Sat, 4 Apr 2015 17:56:01 -0700 Subject: [PATCH 4/8] Removed inaccurate comment. --- geopandas/tools/sjoin.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index 918bc34..eb4f846 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -14,8 +14,6 @@ def sjoin(left_df, right_df, how='left', op='intersects', convert_crs=True, lsuf inner -> use intersection of keys from both dfs; retain only left_df geometry column op: binary predicate {'intersects', 'contains', 'within'} see http://toblerity.org/shapely/manual.html#binary-predicates - use_sindex : Use the spatial index to speed up operation? Default is True - kwargs: passed to op method """ # CHECK VALIDITY OF JOIN TYPE From 2c5dc65786e904349111338658d6923f749028fb Mon Sep 17 00:00:00 2001 From: mdbartos Date: Sat, 4 Apr 2015 18:27:12 -0700 Subject: [PATCH 5/8] Fixed inconsistent indentation --- geopandas/tools/sjoin.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index eb4f846..1511776 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -44,8 +44,6 @@ def sjoin(left_df, right_df, how='left', op='intersects', convert_crs=True, lsuf right_df = right_df.to_crs(left_df.crs) elif left_df.values.nbytes < right_df.values.nbytes: left_df = left_df.to_crs(right_df.crs) - else: - pass # CONSTRUCT SPATIAL INDEX FOR RIGHT DATAFRAME tree_idx = rtree.index.Index() From 54fc5b43f4d1dacd6343346a8a8b4a723099dc07 Mon Sep 17 00:00:00 2001 From: Matthew Bartos Date: Thu, 9 Apr 2015 23:25:26 -0700 Subject: [PATCH 6/8] Removed convert_crs. Added lsuffix and rsuffix to docstring. Compliance with PEP8 and contribution guidelines. Changed default join type to inner. --- geopandas/tools/sjoin.py | 12 ++++-------- tests/test_sjoin.py | 18 +++++++++--------- 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index 1511776..e1a43ad 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -1,10 +1,10 @@ -import geopandas as gpd import numpy as np import pandas as pd import rtree from shapely import prepared +import geopandas as gpd -def sjoin(left_df, right_df, how='left', op='intersects', convert_crs=True, lsuffix='left', rsuffix='right', **kwargs): +def sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuffix='right', **kwargs): """Spatial join of two GeoDataFrames. left_df, right_df are GeoDataFrames @@ -14,6 +14,8 @@ def sjoin(left_df, right_df, how='left', op='intersects', convert_crs=True, lsuf inner -> use intersection of keys from both dfs; retain only left_df geometry column op: binary predicate {'intersects', 'contains', 'within'} see http://toblerity.org/shapely/manual.html#binary-predicates + lsuffix: suffix to apply to overlapping column names (left GeoDataFrame) + rsuffix: suffix to apply to overlapping column names (right GeoDataFrame) """ # CHECK VALIDITY OF JOIN TYPE @@ -38,12 +40,6 @@ def sjoin(left_df, right_df, how='left', op='intersects', convert_crs=True, lsuf # CONVERT CRS IF NOT EQUAL if left_df.crs != right_df.crs: print('Warning: CRS does not match!') - if convert_crs == True: - print('Converting CRS...') - if left_df.values.nbytes >= right_df.values.nbytes: - right_df = right_df.to_crs(left_df.crs) - elif left_df.values.nbytes < right_df.values.nbytes: - left_df = left_df.to_crs(right_df.crs) # CONSTRUCT SPATIAL INDEX FOR RIGHT DATAFRAME tree_idx = rtree.index.Index() diff --git a/tests/test_sjoin.py b/tests/test_sjoin.py index 34c0d4a..f70b280 100644 --- a/tests/test_sjoin.py +++ b/tests/test_sjoin.py @@ -27,7 +27,7 @@ class TestSpatialJoin(unittest.TestCase): shutil.rmtree(self.tempdir) def test_sjoin_left(self): - df = sjoin(self.pointdf, self.polydf, convert_crs=False) + df = sjoin(self.pointdf, self.polydf, how='left') self.assertEquals(df.shape, (21,8)) for i, row in df.iterrows(): self.assertEquals(row.geometry.type, 'Point') @@ -36,8 +36,8 @@ class TestSpatialJoin(unittest.TestCase): def test_sjoin_right(self): # the inverse of left - df = sjoin(self.pointdf, self.polydf, how="right", convert_crs=False) - df2 = sjoin(self.polydf, self.pointdf, how="left", convert_crs=False) + df = sjoin(self.pointdf, self.polydf, how="right") + df2 = sjoin(self.polydf, self.pointdf, how="left") self.assertEquals(df.shape, (12, 8)) self.assertEquals(df.shape, df2.shape) for i, row in df.iterrows(): @@ -46,32 +46,32 @@ class TestSpatialJoin(unittest.TestCase): self.assertEquals(row.geometry.type, 'MultiPolygon') def test_sjoin_inner(self): - df = sjoin(self.pointdf, self.polydf, how="inner", convert_crs=False) + df = sjoin(self.pointdf, self.polydf, how="inner") self.assertEquals(df.shape, (11, 8)) def test_sjoin_op(self): # points within polygons - df = sjoin(self.pointdf, self.polydf, how="left", op="within", convert_crs=False) + df = sjoin(self.pointdf, self.polydf, how="left", op="within") self.assertEquals(df.shape, (21,8)) self.assertAlmostEquals(df.ix[1]['Shape_Leng'], 330454.175933) # points contain polygons? never happens so we should have nulls - df = sjoin(self.pointdf, self.polydf, how="left", op="contains", convert_crs=False) + df = sjoin(self.pointdf, self.polydf, how="left", op="contains") self.assertEquals(df.shape, (21, 8)) self.assertTrue(np.isnan(df.ix[1]['Shape_Area'])) def test_sjoin_bad_op(self): # AttributeError: 'Point' object has no attribute 'spandex' self.assertRaises(ValueError, sjoin, - self.pointdf, self.polydf, how="left", op="spandex", convert_crs=False) + self.pointdf, self.polydf, how="left", op="spandex") def test_sjoin_duplicate_column_name(self): pointdf2 = self.pointdf.rename(columns={'pointattr1': 'Shape_Area'}) - df = sjoin(pointdf2, self.polydf, how="left", convert_crs=False) + df = sjoin(pointdf2, self.polydf, how="left") self.assertTrue('Shape_Area_left' in df.columns) self.assertTrue('Shape_Area_right' in df.columns) @unittest.skip("Not implemented") def test_sjoin_outer(self): - df = sjoin(self.pointdf, self.polydf, how="outer", convert_crs=False) + df = sjoin(self.pointdf, self.polydf, how="outer") self.assertEquals(df.shape, (21,8)) From d43d5e37ef45f21c4ba98592c8441b032d13c439 Mon Sep 17 00:00:00 2001 From: Matthew Bartos Date: Sat, 11 Apr 2015 15:20:21 -0700 Subject: [PATCH 7/8] Better construction of r_idx --- geopandas/tools/sjoin.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index e1a43ad..029ee64 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -52,7 +52,7 @@ def sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuff idxmatch = idxmatch[idxmatch.str.len() > 0] r_idx = np.concatenate(idxmatch.values) - l_idx = np.concatenate((idxmatch.str.len()*pd.Series([[i] for i in idxmatch.index], index=idxmatch.index)).values) + l_idx = idxmatch.index.values.repeat(idxmatch.str.len().values) # VECTORIZE PREDICATE OPERATIONS def find_intersects(a1, a2): From 2af7022ecaff86ea32c3c87eff9fad9345ddaf66 Mon Sep 17 00:00:00 2001 From: Matthew Bartos Date: Sat, 11 Apr 2015 16:07:54 -0700 Subject: [PATCH 8/8] Fixed column widths to conform to PEP8. --- geopandas/tools/sjoin.py | 61 +++++++++++++++++++++++++++++++++------- 1 file changed, 51 insertions(+), 10 deletions(-) diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py index 029ee64..c8d7fe5 100644 --- a/geopandas/tools/sjoin.py +++ b/geopandas/tools/sjoin.py @@ -4,14 +4,16 @@ import rtree from shapely import prepared import geopandas as gpd -def sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuffix='right', **kwargs): +def sjoin(left_df, right_df, how='inner', op='intersects', + lsuffix='left', rsuffix='right', **kwargs): """Spatial join of two GeoDataFrames. left_df, right_df are GeoDataFrames how: 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 + inner -> use intersection of keys from both dfs; + retain only left_df geometry column op: binary predicate {'intersects', 'contains', 'within'} see http://toblerity.org/shapely/manual.html#binary-predicates lsuffix: suffix to apply to overlapping column names (left GeoDataFrame) @@ -48,7 +50,8 @@ def sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuff tree_idx.insert(i, right_df_bounds[i]) # FIND INTERSECTION OF SPATIAL INDEX - idxmatch = left_df['geometry'].apply(lambda x: x.bounds).apply(lambda x: list(tree_idx.intersection(x))) + idxmatch = (left_df['geometry'].apply(lambda x: x.bounds) + .apply(lambda x: list(tree_idx.intersection(x)))) idxmatch = idxmatch[idxmatch.str.len() > 0] r_idx = np.concatenate(idxmatch.values) @@ -61,27 +64,65 @@ def sjoin(left_df, right_df, how='inner', op='intersects', lsuffix='left', rsuff def find_contains(a1, a2): return a1.contains(a2) - predicate_d = {'intersects': find_intersects, 'contains': find_contains, 'within': find_contains} + predicate_d = {'intersects': find_intersects, + 'contains': find_contains, + 'within': find_contains} check_predicates = np.vectorize(predicate_d[op]) # CHECK PREDICATES - result = pd.DataFrame(np.column_stack([l_idx, r_idx, check_predicates(left_df['geometry'].apply(lambda x: prepared.prep(x)).values[l_idx], right_df['geometry'].values[r_idx])])) + result = ( + pd.DataFrame( + np.column_stack( + [l_idx, + r_idx, + check_predicates( + left_df['geometry'] + .apply(lambda x: prepared.prep(x)).values[l_idx], + right_df['geometry'].values[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) + result = ( + pd.DataFrame(result[result['match_bool']==1]) + .drop('match_bool', axis=1) + ) # IF 'WITHIN', SWAP NAMES AGAIN 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)}) + result = result.rename(columns={ + 'index_%s' % (lsuffix): 'index_%s' % (rsuffix), + 'index_%s' % (rsuffix): 'index_%s' % (lsuffix)}) # APPLY JOIN 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('geometry', axis=1), left_on='index_%s' % rsuffix, right_index=True, suffixes=('_%s' % lsuffix, '_%s' % rsuffix)) + return ( + left_df + .merge(result, left_index=True, right_index=True) + .merge(right_df.drop('geometry', 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('geometry', axis=1), how='left', left_on='index_%s' % rsuffix, right_index=True, suffixes=('_%s' % lsuffix, '_%s' % rsuffix)) + return ( + left_df + .merge(result, left_index=True, right_index=True, how='left') + .merge(right_df.drop('geometry', axis=1), + how='left', left_on='index_%s' % rsuffix, right_index=True, + suffixes=('_%s' % lsuffix, '_%s' % rsuffix)) + ) elif how == 'right': - return left_df.drop('geometry', 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) + return ( + left_df + .drop('geometry', 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) + )