3

Coworker and I are getting different results from the same code on the same data. It is a very simple script where the shapefile is imported and converted to 4326, then we import the points and do a spatial join.

import pandas as pd
import geopandas as gpd

shp = gpd.read_file("shapefile.shp")
shp = shp.to_crs("EPSG:4326")

points = pd.read_csv("points.csv")
gdf = gpd.GeoDataFrame(points, crs="EPSG:4326", geometry=gpd.points_from_xy(points.lon, points.lat))

sjoin = gdf.sjoin(shp, how='inner')

When I run it, I get 59/736 points returned, but my coworker only gets 3/736.

We are both reading in the same points CSV and shapefile from the server and running the exact same code. We restarted our computers, created new environments, and confirmed our versions of GeoPandas, rtree, etc. are all the same. But we still get different results.

Has anyone run into this or something similar? Not sure what to do about it.

-- Update: Environments -- My coworker is off today, so I can't get a full 'conda list' from them, but we did previously confirm we have the same versions of the following:

  • python 3.10.1 h9a09f29_2_cpython conda-forge
  • geopandas 0.10.2 pyhd8ed1ab_0 conda-forge
  • geopandas-base 0.10.2 pyha770c72_0 conda-forge
  • pyproj 3.3.0 py38h12f9802_0 conda-forge
  • rtree 0.9.7 py310h1cbd46b_3 conda-forge

We built our environments following the exact steps below:

conda create --name env2
conda install -c conda-forge geopandas

Something interesting though was that the above commands originally returned slightly different environments. My environment was built with pyproj 3.3.0 and their environment was built with pyproj 2.6.1.post1. We thought this could be the problem, but updating to 3.3.0 did not change their results.

-- Update: Data -- We both get the same output for points.shape, 736x8. I can't currently confirm we get the same output for shp.shape, but the count of polygons is very important and they would have raised a flag if it was different from what it should always be. Though I will update once they get back to me.

The 3 points they get are not points I get in my 59.

    geometry
284 POINT (-155.47659 19.20657)
285 POINT (-157.72620 21.28898)
286 POINT (-157.72621 21.28953)

For some context, these are all points falling very close to the border of a polygon. On my end, these points are not considered to be intersecting the polygon, but on coworker's end they were.

enter image description here

10
  • shp.shape and points.shape give the same output for the both of you?
    – BERA
    Dec 30, 2021 at 18:54
  • @BERA yes, we get the same shapes
    – chchchchch
    Dec 30, 2021 at 19:42
  • 1
    Can you share your data with us? Can you please show the result of print(gdf) just before the join procedure?
    – Taras
    Dec 30, 2021 at 21:33
  • 2
    can you also provide the output of pip freeze and python -V for both computers ? even if both are the same it would be helpful for triaging exactly which version of the code ks being run
    – jaksco
    Dec 31, 2021 at 5:23
  • @Taras sorry, I am not able to share the data at the moment, but I can check.
    – chchchchch
    Dec 31, 2021 at 16:46

1 Answer 1

1

After a few days of confusion, we figured it out it was probably being caused by a caching issue between Jupyter Lab and Jupyter Notebook. Once the notebook had been run in Jupyter Lab and returned 59 points, creating a new environment or running in Jupyter Notebook would not change the output. However, if I created a new environment without defaults, only installed geopandas and then installed Notebook through Anaconda Navigator, I would get the 3 points.

At the end of the day I still don't know the exact cause, by I felt the easiest fix was to just start from scratch. I uninstalled Anaconda and Python, and re-installed. Now we are getting the same results.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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