Projection conversion with Python/Pyproj: Canada Albers Equal Area to WGS

Question

How do I get the geographical coordinates (lat and lon) from x,y values in Canada Albers Equal Area projection?

The range of x and y values as follows:

min(x) = -5171461.906673405 max(x) = 3372679.8809247203 min(y) = 209275.49942131247 max(y) = 5395930.584568907

I tried specifying EPSG code of Canada Equal Area as 102001 (see https://epsg.io/102001), but this throws an error in Pyproj.

RuntimeError: no options found in 'init' file

Gdalinfo for the tiff files of projected coordinate system and WGS are provided below. If you see the upper right corner value (0d38'48.11"W) in geographic coordinates of the projected coordinate system gdalinfo output, it looks like something is wrong -- totally off.

Coordinate System is:
PROJCS["Canada_Albers_Equal_Area_Conic",
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.2572221010042,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4269"]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["standard_parallel_1",50],
PARAMETER["standard_parallel_2",70],
PARAMETER["latitude_of_center",40],
PARAMETER["longitude_of_center",-96],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
Origin = (-5171461.906673405319452,5395930.584568906575441)
Pixel Size = (12034.002517743838325,-12034.002517743838325)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-5171461.907, 5395930.585) (165d 9'48.70"E, 46d33'32.99"N)
Lower Left  (-5171461.907,  197241.497) (145d24'33.39"W, 22d41' 7.36"N)
Upper Right ( 3384713.883, 5395930.585) (  0d38'48.11"W, 62d22'40.07"N)
Lower Right ( 3384713.883,  197241.497) ( 60d 4'49.15"W, 33d 7'48.19"N)
Center      ( -893374.012, 2796586.041) (114d45'33.46"W, 63d59'36.53"N)
Band 1 Block=711x2 Type=Float32, ColorInterp=Gray
NoData Value=-999

Gdalinfo of the WGS file is here

Driver: GTiff/GeoTIFF
Files: ABoVE.hydro.wy.OUT_PREC.1987.tif
Size is 476, 166
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-172.125000000000000,41.875000000000000)
Pixel Size = (0.249474789915966,0.248493975903614)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-172.1250000,  41.8750000) (172d 7'30.00"W, 41d52'30.00"N)
Lower Left  (-172.1250000,  83.1250000) (172d 7'30.00"W, 83d 7'30.00"N)
Upper Right ( -53.3750000,  41.8750000) ( 53d22'30.00"W, 41d52'30.00"N)
Lower Right ( -53.3750000,  83.1250000) ( 53d22'30.00"W, 83d 7'30.00"N)
Center      (-112.7500000,  62.5000000) (112d45' 0.00"W, 62d30' 0.00"N)
Band 1 Block=476x4 Type=Float32, ColorInterp=Gray
NoData Value=-9999

Answer 1

Canada Albers Equal Area Conic: ESRI:102001 is not an EPSG SRID value but an ESRI SRID value not recognized by PyProj.

So as FSimardGIS says, simply use the PROJ4 string of Canada Albers Equal Area Conic: ESRI:102001.

import pyproj
# ESRI:102001 proj4 string
inproj = pyproj.Proj('+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs')
# Wgs84
outproj = pyproj.Proj(init='epsg:4326')
xmin,ymin = (5171461.906673405,209275.49942131247)
pyproj.transform(inproj,outproj,xmin,ymin)
(-46.520268732851214, 22.775277555142978)
xmax,ymax = (3372679.8809247203,5395930.584568907)
pyproj.transform(inproj,outproj,xmax,ymax)
(-0.682414393773495, 62.48292643585552)

If you want to get the PROJ4 string automatically, you can use PyCRS:

import pycrs
crs = pycrs.parse.from_esri_code(102001)
inproj = pyproj.Proj( crs.to_proj4())
pyproj.transform(inproj,outproj,xmin,ymin)
(-46.520268732851214, 22.775277555142978)

Answer 2

pyproj 2.2.0 can solve this:

>>> import pyproj
>>> pyproj.__version__
'2.2.0'
>>> transformer = pyproj.Transformer.from_crs("esri:102001", "epsg:4326")
>>> transformer.transform((-5171461.906673405, 3372679.8809247203), (209275.49942131247, 5395930.584568907))
((22.775277555142978, 62.48292643585552), (-145.4797312671488, -0.682414393773495))