jGaboardi / tigernet

Network topology via TIGER/Line Edges
https://jgaboardi.github.io/tigernet/
BSD 3-Clause "New" or "Revised" License
6 stars 1 forks source link

new pyproj causes variations in projection results #67

Closed jGaboardi closed 3 years ago

jGaboardi commented 3 years ago

Something changed between pyproj 2.6.1.post1 and pyproj 3.0.0.post1 that is leading to small variations in projected coordinates, which is then propagated on to point snapping distance, cost matrices, etc. For now, pinning to 2.6.1.post1 is reasonable (#66). A longer term solution is either updating the tests to reflect the altered projection results or figuring out if this is an actual bug. The code chunk below demonstrates the output difference between pyproj 2.6.1.post1 (top) and pyproj 3.0.0.post1 (bottom) with an example through geopandas and a MWE with only pyproj.

@knaaptime @slumnitz @martinfleis Have any of you noticed this in your work?

import warnings
warnings.filterwarnings("ignore")
import geopandas
import numpy
import pyproj
import shapely

packages = [geopandas, numpy, pyproj, shapely]
for p in packages:
    print(f"{p.__name__}: {p.__version__}")

layer = "WeightedParcels_Leon_FL_2010"
fin = f"zip://{layer}.zip!{layer}.shp"
bbox = (-84.279, 30.480, -84.245, 30.505)
gdf = geopandas.read_file(fin, bbox=bbox)
TEST_PARCEL = "1117160000020"
gdf = gdf[gdf["PARCEL_ID"] == TEST_PARCEL].reset_index(drop=True)
point = gdf.geometry[0]
x, y = point.x, point.y

from_crs, to_crs = gdf.crs, "EPSG:2779"

print(f" -- {geopandas.__name__}/{shapely.__name__} --")
print(gdf.crs, point)
gdf = gdf.to_crs(to_crs)
point = gdf.geometry[0]
print(gdf.crs, point)

print(f" -- {pyproj.__name__} --")
transformer = pyproj.Transformer.from_crs(from_crs, to_crs)
print(from_crs, (x, y))
print(to_crs, transformer.transform(y, x))
geopandas: 0.8.1
numpy: 1.19.4
pyproj: 2.6.1.post1
shapely: 1.7.1
 -- geopandas/shapely --
epsg:4326 POINT (-84.25873245332716 30.48436516338524)
EPSG:2779 POINT (623164.270749338 164564.2000569425)
 -- pyproj --
epsg:4326 (-84.25873245332716, 30.48436516338524)
EPSG:2779 (623164.270749338, 164564.20005694253)
geopandas: 0.8.2
numpy: 1.20.1
pyproj: 3.0.0.post1
shapely: 1.7.1
 -- geopandas/shapely --
epsg:4326 POINT (-84.25873245332716 30.48436516338524)
EPSG:2779 POINT (623164.6468275142 164563.5686168194)
 -- pyproj --
epsg:4326 (-84.25873245332716, 30.48436516338524)
EPSG:2779 (623164.6468275142, 164563.56861681942)
jGaboardi commented 3 years ago

See original issue in #65

martinfleis commented 3 years ago

I think that this is a result of changes in how proj handles transformation grids. The new result should have the higher precision compared to the other one, so I would fix the tests. There was one issue of this kind in geopandas some time ago.