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).
New info, code works with gdal python library 2.2.2 but not 3.0.2
"""
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
