3

I have a GeoTIFF file and I wanted to know its CRS. When I open it on QGIS I get:GeoTIFF descriptions.

I wanted to see if I could get the same result on Python (EPSG = 5936, the code for Alaska_Polar_Stereographic), I used:

ds = ds = gdal.Open('/content/drive/MyDrive/MB_Model/Glacier_Mask.tif')
proj = osr.SpatialReference(wkt=ds.GetProjection())
print(proj.GetAttrValue('AUTHORITY',1))

and got as answer: 4326

I also used: !gdalinfo '/content/drive/MyDrive/MB_Model/Glacier_Mask.tif'

and got:

Driver: GTiff/GeoTIFF
Files: /content/drive/MyDrive/MB_Model/Glacier_Mask.tif
Size is 523, 582
Coordinate System is:
PROJCS["WGS_1984_EPSG_Alaska_Polar_Stereographic",
    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"]],
    PROJECTION["Stereographic"],
    PARAMETER["latitude_of_origin",90],
    PARAMETER["central_meridian",-150],
    PARAMETER["scale_factor",0.994],
    PARAMETER["false_easting",2000000],
    PARAMETER["false_northing",2000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (2522546.904911093413830,-1317885.149017794057727)
Pixel Size = (125.524826080368911,-125.524826080368996)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 2522546.905,-1317885.149) (141d 2'59.20"W, 60d23'49.86"N)
Lower Left  ( 2522546.905,-1390940.598) (141d14'22.49"W, 59d47'19.46"N)
Upper Right ( 2588196.389,-1317885.149) (139d56'49.28"W, 60d18'20.62"N)
Lower Right ( 2588196.389,-1390940.598) (140d 9'33.58"W, 59d41'58.00"N)
Center      ( 2555371.647,-1354412.873) (140d35'56.88"W, 60d 2'56.06"N)
Band 1 Block=523x3 Type=Float32, ColorInterp=Gray
  NoData Value=-9999
  Metadata:
    SourceBandIndex=0

I am completely lost because the three different methods are giving me three different results. Why isn't the code 5936 appearing anywhere in the gdalinfo command ?

5
  • I don't see 5901 in any of your results - it looks like gdalinfo and qgis agree, not sure whats wrong with your python code though
    – Ian Turton
    Nov 5 at 20:10
  • 1
    Perhaps print(proj.GetAttrValue('AUTHORITY',1)) prints the authority code of the datum DATUM["WGS_1984" as it appears in the gdalinfo report.
    – user30184
    Nov 5 at 21:06
  • Sorry my mistake, the EPSG for Alaska Polar Stereographic is 5936. I'm just surprised it doesn't appear on the gdalinfo, it just gives the projection without its code.
    – vdc
    Nov 5 at 22:35
  • 2
    QGIS isn't giving you the EPSG code either based on that screenshot. The projection WKT stored in the GeoTIFF metadata doesn't include a PROJCS AUTHORITY, so GDAL falls back to giving you the GEOGCS AUTHORITY.
    – user2856
    Nov 6 at 5:13
  • Ok so that might be why then, thanks for the insight !
    – vdc
    Nov 7 at 23:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.