Unable to attach a projection before creating a GDAL image

Question

I'm unable to add a projection to a newly-created raster with GDAL-python, using the code below:

arr = np.random.rand(600, 400)
filename = 'path_to_file_'
driver = gdal.GetDriverByName('GTiff')
x_size = arr.shape[1]
y_size = arr.shape[0]

dataset = driver.Create(filename, x_size, y_size,
                        eType=gdal.GDT_Float32)
dataset.SetProjection(projection)
dataset.GetRasterBand(1).WriteArray(arr)

It's expected in the code above that the projection isn't attached to the dataset, because the dataset has already been written when the projection is added to it. Is there another way to create an image in a way that the projection is attached to the saved image?

Note that the projection (WGS84) that I'm trying to attach to the image to create was extracted from an image used as input to create the newly-created image (i.e. I'm just trying to keep the projection of the source image in the image to save). Here is the string returned by dataset.GetProjection() (i.e. the content of the variable projection):

'GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","4326"]]'

Answer 1

Your code works fine. Make sure you flush the writes to disk by dereferencing the dataset variable (a well known python GDAL "gotcha").

Working example:

from osgeo import gdal
import numpy as np

arr = np.random.rand(600, 400)
filename = 'path_to_file.tif'
driver = gdal.GetDriverByName('GTiff')
x_size = arr.shape[1]
y_size = arr.shape[0]
projection = 'GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwicbashh",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433],AXIS["Geodetic longitude",EAST],AXIS["Geodetic latitude",NORTH],AUTHORITY["EPSG","4326"]]'

dataset = driver.Create(filename, x_size, y_size,
                        eType=gdal.GDT_Float32)
dataset.SetProjection(projection)
dataset.GetRasterBand(1).WriteArray(arr)

del dataset # dereference the variable so writes get flushed to disk 

gdalinfo output:

gdalinfo path_to_file.tif 
Driver: GTiff/GeoTIFF
Files: path_to_file.tif
Size is 400, 600
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"]]
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  600.0)
Upper Right (  400.0,    0.0)
Lower Right (  400.0,  600.0)
Center      (  200.0,  300.0)
Band 1 Block=400x5 Type=Float32, ColorInterp=Gray