0

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
Improve this question
1
  • 1
    Can you post the code that you used to perform the transformation? Perhaps using the projection string could be a way to solve the issue, like this : srcProj = "+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m no_defs"
    – FSimardGIS
    Jan 5, 2019 at 2:38

2 Answers 2

Reset to default
2

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)
Improve this answer
1
  • Thanks for the quick answer! I tried this code. It seems to work except that it has one typo (pycrs.parse"r") and in the next line 'inproj=pyproj.Proj(crs.to_proj4())', I get an error 'RuntimeError: conic lat_1 = -lat_2', but the problem seems to be something else. I've added more info to the question detailing it, could you please have a look and let me know if you spot anything?
    – Solomon Vimal
    Jan 8, 2019 at 0:36
1

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))
Improve this answer

Your Answer

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

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