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? Commented Apr 18, 2017 at 20:28
  • @hakim you can ask another question. Commented 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? Commented 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. Commented 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. Commented May 3, 2017 at 9:28

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.