0

I'm trying to compute a CHM (Canopy Height Model) using GDAL for Python.

A_TIF = os.getcwd() + xyzlocation + folder + "\\dtm.tif"
B_TIF = os.getcwd() + xyzlocation + folder + "\\dsm.tif"
OUT_TIF = os.getcwd() + xyzlocation + folder + "\\chm.tif"

ds1 = gdal.Open(A_TIF, gdalconst.GA_ReadOnly)
ds2 = gdal.Open(B_TIF, gdalconst.GA_ReadOnly)

arr1 = ds1.ReadAsArray()
arr2 = ds2.ReadAsArray()

# subtract every value from both rasters
result = np.subtract(arr2, arr1)

# get attributes from input raster 1 and set them to the output raster
driver = gdal.GetDriverByName("GTiff")
output = driver.Create(OUT_TIF, arr1.shape[1], arr1.shape[0], 1, gdal.GDT_Float32)
output.SetGeoTransform(ds1.GetGeoTransform())
output.SetProjection("EPSG:3067") # (Error) GDAL failure (1) missing [

output.GetRasterBand(1).WriteArray(result)

The line mentioned above produces an error:

(Error) GDAL failure (1) missing [

Why is this?


Also, this code may fail to produce the CHM (results in 1KB file). But this is not the point of the question.

2

Try:

from osgeo import osr
proj = osr.SpatialReference()
proj.ImportFromEPSG(3067)
output.SetProjection(proj)
|improve this answer|||||
2

Building on @Yjelza's answer, the SetProjection() method takes a string as its argument and not a SpatialReference object. You have to call the ExportToWkt() method on the object.

from osgeo import osr
sr = osr.SpatialReference()
sr.ImportFromEPSG(3067)
output.SetProjection(sr.ExportToWkt())

If you want to take the projection of one of the input rasters it would be as simple as you did by setting the GeoTransform.

output.SetProjection(ds1.GetProjection())

Note that if you are trying to use a different projection that the one ds1 has, the GeoTransform object won't work because it has the raster origin coordinates in another spatial reference.


Also, it is possible that running your script results in a 1KB file because you are not closing the file at the end of the script.

Try adding the following line at the end:

output = None
|improve this answer|||||
  • 1
    Good spot, I forgot to ExportToWkt! But I would suggest output = None instead of del, it helps with properly closing the dataset in the gdal bindings – Yjelza Feb 26 at 16:10
  • @Yjelza Didn't know there was a difference between using del and setting it to None. Will keep that in mind! – Marcelo Villa-Piñeros Feb 26 at 16:19

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.