2

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"]]'

Improve this question
0

1 Answer 1

Reset to default
2

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
Improve this answer
7
  • It's probably working but I can't verify it for the time being, because the image I'm generating with SNAP seems to lack the right ellipsoid (even when I'm terrain-correcting it with WGS84 as the map projection). I noticed this problem both with projection.GetProjection and with QGIS. Any idea what's causing this?
    – Hakim
    Apr 18, 2017 at 20:28
  • @hakim you can ask another question.
    – user2856
    Apr 18, 2017 at 23:03
  • I've tested your code with and without del dataset, and the projection was the same both with gdalinfo and with dataset.GetProjection (I'm using gdal 2.1.2). I think SNAP is unable to read the coordinates because of the Corner Coordinates parameter, which is expressed in pixels instead of being given as a (lat, lon) angle in degrees. Do you know how to convert this parameter to degrees?
    – Hakim
    Apr 26, 2017 at 22:26
  • @hakim Using del to dereference the dataset is only required if you need to access the dataset before your script finishes (or you use an IDE that holds the variables open after the script finishes). I don't know why snap doesn't understand the gdal georeferencing.
    – user2856
    Apr 27, 2017 at 7:14
  • 1
    @hakim yes that is correct. If you had asked about attaching georeferencing as well as projection info I would have included that in the example.
    – user2856
    May 3, 2017 at 9:28

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.