1

I have two geodataframes that I want to perform a spatial join on; one is points, one is polygons. They are both dtype: int64. They both have .crs = {'init': 'epsg:4326'}. I have verified both of them in ArcMap that the points are indeed within the polys. I even performed an intersect between these points and polygons in ArcMap and it produced the correct new feature class (I have reasons for wanting to do this without arpy/manually in ArcMap). **Note: the points are from a shapefile and the polys are from a feature class enter image description here

However, when I turn my centroid points and parcel polygons into geodataframes and use .sjoin(), it returns an empty geodataframe. I want the output to be a point-based geometry so that I can eventually turn it back into a point shapefile. I have tried:

intersect_test = gpd.sjoin(cents, polys, how='left', op='intersects')

and

intersect_test = gpd.sjoin(polys, cents, how='inner', op='intersects')

and

intersect_test = gpd.sjoin(polys, cents, how='right', op='intersects')

and

intersect_test = gpd.sjoin(cents, polys, how='left', op='intersects')

and pretty much every other configuration I can think of and it either returns the data for one of them, or returns a completely empty geodataframe result.

 intersect_test = gpd.sjoin(cents, polys, how='inner', op='intersects')
 intersect_test.count()
Out[126]: 
BLD_UNITS       0
LAND_USE_T      0
PARCEL_ID       0
PROP_IND_T      0
STORY_NBR       0
geometry        0
index_right     0
GEOID           0
CensusPop       0
CBArea          0
ST_FIPS         0
Shape_Length    0
Shape_Area      0
CO_FIPS         0
HU_Pop          0
Sq_Ft           0
dtype: int64

How can I resolve this without having to manually perform the intersect in ArcMap?

EDIT Here's how I made my geodataframes. As I mentioned; one is a shapefile and one is a .gdb feature class. Question: Could this be due to crs-related problem?

#create poly gdf
cb_gdb = r"C:\Projects\Pop_Alloc\CensusBlocksStates.gdb"
cb_feat = "CBs_{}".format(state)
cents = gpd.read_file(cb_gdb, layer=cb_feat)
cents = {'init': 'epsg:4326'}

#create point gdf. The parcel centroids are first created from a polygon 
gdf. The new gdf is written to a shapefile to be tested against the census 
block polygons (to make sure they do in fact fall within the boundaries of 
the census blocks (cbs)) and the new centroid gdf(its a 
`type='geopandas.geodataframe.GeoDataFrame'`) is then used in the `.sjoin`       
cents = parcel_res_df
cents['geometry'] = 
cents['geometry'].representative_point()
cents_out_feat = r"C:\Projects\Pop_Alloc\{}_Res_centroids.shp".format(state)
cents = {'init': 'epsg:4326'}
cents.to_file(cents_out_feat)
  • 1
    What you are trying should work. Can you show how you create the dataframes? – BERA Jun 27 at 20:18
  • Just edited to include how I made the .gdf's. – gwydion93 Jun 28 at 13:23
  • I would try plotting them in geopandas and see if they overlap. – BERA Jun 28 at 13:45
  • try it with overlay from geopandas geopandas.org/set_operations.html – Leo Jun 28 at 13:50
  • Tried overlay but got: TypeError: overlay only takes GeoDataFrames with (multi)polygon geometries. I have polys and points. – gwydion93 Jun 28 at 14:43
0

This problem is a coordinate reference system issue. It is most likely caused by my census block feature class(aka polys) being in ESRI:102008. When I converted my polys feature class to a geodataframe, it showed {} when I ran polys.crs. Further investigation showed that the geometry looked like : (POLYGON ((1385968.793600001 83806.20800000057..., which was also incorrect. The solution required first finding the original crs, assigning the 'polys.crs' to the correct crs and then reprojecting to the desired crs. To do this, I used ogr/gdal:

import ogr
def get_proj(gdb_input, fc ):
    driver = ogr.GetDriverByName("OpenFileGDB")
    gdb = driver.Open(gdb_input, 0)
    print('There are {} total layers in {}'.format(gdb.GetLayerCount(), gdb_input))

    for i in range(gdb.GetLayerCount()):
        layer = gdb.GetLayerByIndex(i)
        if layer.GetDescription() == fc:
            break
    in_crs = layer.GetSpatialRef().ExportToProj4()
    return in_crs
    print('CRS = ', in_crs)

polys.crs = get_proj(poly_gdb, poly_name)
polys = cb_df.to_crs(cents.crs)
intersect = gpd.sjoin(ccents, polys, how='inner', op='intersects')

This returned the correctly intersected geodataframe. I have to say that I am surprised there is not documentation for what seems like a very common problem in using Geopandas with data from an ESRI context. Hope this helps.

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.