Map projections transform of point feature using different versions of Pyproj

In many cases, we may need to change the projection of geographic features. This is a very common process in GIS and other spatial analyses. You can use QGIS or ArcGIS to change the projection of shapefile pretty easily. However, in some cases, you may need to change the project of point features using programming. Pyproj is a very popular tool for changing the projection of different point features. However, the different versions of Pyproj have very different call functions to do the projection transform.

Here is a sample code I wrote to do the projection transform using Python with consideration of different versions of Pyproj.

def reprojPnt(lon, lat, epsg_code):
    '''This function is used to reproject the lon, lat
    parameters:
        lon: the longituted, float
        lat: the latitude, epsg 4326, float
        epsg_code: the output epsg code, int
    return:
        (reproj_lon, reproj_lat)
    
    @Xiaojiang Li, Temple University, August 22, 2020
    '''
    
    import pyproj
    from shapely.geometry import Point, mapping, shape
    from functools import partial
    from shapely.ops import transform
    
    
    point_degree = Point(float(lon), float(lat))
    
    # for different version of pyproject
    if pyproj.__version__[0]=='1': #proj1.9.6, you may have error to reproject lon, lat to mass proj 6492
        project = partial(pyproj.transform,pyproj.Proj(init='EPSG:4326'),pyproj.Proj(init='EPSG:'+str(epsg_code))) #3857 is psudo WGS84 the unit is meter
        point_meter = transform(project, point_degree)
    else:
        project = pyproj.Transformer.from_crs(4326, epsg_code, always_xy=True)
        x2, y2 = project.transform(lon, lat)
        # print('The new coordinate is:', x2, y2)
        point_meter = Point(float(x2), float(y2))
        # print('the point is:', point_meter)
    
    return point_meter.coords[0]

Then you can call the above function to reproject the point feature,

#6492 is the epsg code for local projection of Massachusetts
reprojPnt(-71.12704, 42.387634, 6492)

Here is a tricky thing. If you are using Pyproj-1.9.6, you will get coordinates,

(230710.7523346356, 904192.2213170032)

It looks right. However, if you double-check the street map projected in EPSG: 6492, which is located in Cambridge, MA, the result is wrong!!! Check the bottom of the figure.

Figure1. The street map in the Boston metro area, MA.

I was struggling with this for a while because it works for EPSG: 3857, which is the WGS 84 / Pseudo-Mercator with the unit of meter. I actually use this project in creating sample points along the street in my Treepedia project.

Then I decided to try with a newer version of Pyproj, hope it would work. I upgraded my pyproj to 2.6.1 by typing in terminal or Anaconda prompt,

pip install pyproj==2.6.1

Then your old pyproj will be uninstalled and install the pyproj-2.6.1. I run the statement again,

reprojPnt(-71.12704, 42.387634, 6492)

Here is the new result,

(756923.526617883, 2966503.97943754)

It works. I suggest upgrading your Pyproj if you don’t want to have this error in the projection transform.