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, arr1.shape, 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.