0

I can't figure out how to get the projection info from a raster file using gdal. I thought my raster had a bad srs so I used a script from # https://gdal.org/tutorials/raster_api_tut.html to create a raster with a correct srs built with gdal. Gdalinfo.py reports a correct srs. This is my code... projection is alway and empty (that should indicate an absent or invalid srs (but this is not the case).

Edit: New info, code works with gdal python library 2.2.2 but not 3.0.2

Thanks for you help!

"""
Created on Sun Oct 25 15:09:05 2020

@author: Nicolas
"""
#https://stackoverflow.com/questions/37648439/simplest-way-to-save-array-into-raster-file-in-python

from osgeo import gdal, osr
import numpy

# https://gdal.org/tutorials/raster_api_tut.html
fileformat = "GTiff"
driver = gdal.GetDriverByName(fileformat)
metadata = driver.GetMetadata()
if metadata.get(gdal.DCAP_CREATE) == "YES":
    print("Driver {} supports Create() method.".format(fileformat))

if metadata.get(gdal.DCAP_CREATECOPY) == "YES":
    print("Driver {} supports CreateCopy() method.".format(fileformat))
dst_ds = driver.Create(r"c:\temp\gdal.tif", xsize=512, ysize=512,
                    bands=1, eType=gdal.GDT_Byte)
dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])
srs = osr.SpatialReference()
srs.SetUTM(11, 1)
srs.SetWellKnownGeogCS("NAD27")
dst_ds.SetProjection(srs.ExportToWkt())
print ('srs = ',srs)# this is good
raster = numpy.zeros((512, 512), dtype=numpy.uint8)
dst_ds.GetRasterBand(1).WriteArray(raster)
# Once we're done, close properly the dataset
srs = None
dst_ds = None #srs is file and well georeferenced in Qgis.

# ---------------------------------------------------------------------------
# open dataset try to read srs
# ---------------------------------------------------------------------------
raster_ds = gdal.Open(r"C:\temp\gdal.tif")
# first try
print("Projection is {}".format(raster_ds.GetProjection()))
# second try
srs = osr.SpatialReference()
srs.ImportFromWkt(raster_ds.GetProjectionRef())
print ('srs =', srs)
# thrid try
srs = osr.SpatialReference()
srs = raster_ds.GetProjection()
print('srs =', srs)
srs = raster_ds.GetProjectionRef()
# forth try
# from gdal_info https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/gdalinfo.py
pszProjection = raster_ds.GetProjectionRef()
print(pszProjection)
if pszProjection is not None:
    hSRS = osr.SpatialReference()
    if hSRS.ImportFromWkt(pszProjection) == gdal.CE_None:
        pszPrettyWkt = hSRS.ExportToPrettyWkt(False)
        print("Coordinate System is:\n%s" % pszPrettyWkt)
    else:
        print("Coordinate System is `%s'" % pszProjection)

# ----------------------------------------------------------------------
# Results for gdal info
# ----------------------------------------------------------------------
# Driver: GTiff/GeoTIFF
# Files: C:/temp/gdal.tif
# Size is 512, 512
# Coordinate System is:
# PROJCRS["NAD27 / UTM zone 11N",
# BASEGEOGCRS["NAD27",
# DATUM["North American Datum 1927",
# ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
# LENGTHUNIT["metre",1]]],
# PRIMEM["Greenwich",0,
# ANGLEUNIT["degree",0.0174532925199433]],
# ID["EPSG",4267]],
# CONVERSION["UTM zone 11N",
# METHOD["Transverse Mercator",
# ID["EPSG",9807]],
# PARAMETER["Latitude of natural origin",0,
# ANGLEUNIT["degree",0.0174532925199433],
# ID["EPSG",8801]],
# PARAMETER["Longitude of natural origin",-117,
# ANGLEUNIT["degree",0.0174532925199433],
# ID["EPSG",8802]],
# PARAMETER["Scale factor at natural origin",0.9996,
# SCALEUNIT["unity",1],
# ID["EPSG",8805]],
# PARAMETER["False easting",500000,
# LENGTHUNIT["metre",1],
# ID["EPSG",8806]],
# PARAMETER["False northing",0,
# LENGTHUNIT["metre",1],
# ID["EPSG",8807]]],
# CS[Cartesian,2],
# AXIS["(E)",east,
# ORDER[1],
# LENGTHUNIT["metre",1]],
# AXIS["(N)",north,
# ORDER[2],
# LENGTHUNIT["metre",1]],
# USAGE[
# SCOPE["unknown"],
# AREA["North America - 120°W to 114°W and NAD27 by country - onshore"],
# BBOX[26.93,-120,78.13,-114]],
# ID["EPSG",26711]]
# Data axis to CRS axis mapping: 1,2
# Origin = (444720.000000000000000,3751320.000000000000000)
# Pixel Size = (30.000000000000000,-30.000000000000000)
# Metadata:
# AREA_OR_POINT=Area
# Image Structure Metadata:
# INTERLEAVE=BAND
# Corner Coordinates:
# Upper Left ( 444720.000, 3751320.000) (117d35'52.47"W, 33d54' 9.25"N)
# Lower Left ( 444720.000, 3735960.000) (117d35'49.00"W, 33d45'50.54"N)
# Upper Right ( 460080.000, 3751320.000) (117d25'54.40"W, 33d54'11.75"N)
# Lower Right ( 460080.000, 3735960.000) (117d25'51.90"W, 33d45'53.03"N)
# Center ( 452400.000, 3743640.000) (117d30'51.94"W, 33d50' 1.25"N)
# Band 1 Block=512x16 Type=Byte, ColorInterp=Gray
1

With the help of the gdal-dev@lists dot osgeo dot org, I was able to figure that this problem was caused by my Anaconda environment. Problem was solved by updating gdal 3.0.2 (the highest version I could install using the conda navigator) to version 3.1.4 using conda install -c conda-forge gdal. Thanks P. H. Nicolas

| improve this answer | |

Your Answer

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

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