Search code examples
pythonnumpygdal

gdal WriteArray() crashes python without a stack trace


I'm trying to write an array to a geotiff using gdal. Each row of the array is identical, and I used np.broadcast_to to create the array.

When I try to write it, the I get a windows popup saying "Python has stopped working: A problem cause the program to stop working correctly. Please close the program"

This approximates the steps I'm taking:

import gdal
import numpy as np

driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create("C:/raster.tif", 1000, 1000, 1, 6)
band = outRaster.GetRasterBand(1)

# Create  array 
a = np.arange(0,1000, dtype='float32')
a1 = np.broadcast_to(a, (1000,1000))

# try writing
band.WriteArray(a1) # crash

Solution

  • The problem is that the input array created by broadcast_to isn't contiguous on disk. As described in the numpy documentation, more than one element array may point to the same memory address. This causes problems in gdal.

    Instead of using broadcast_to, use something that stores each element as its own place in memory. As an illustrative example, see the following code:

    import gdal
    import numpy as np
    import sys
    
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create("C:/raster.tif", 1000, 1000, 1, 6)
    band = outRaster.GetRasterBand(1)
    
    # Create 1000 x 1000 array two different ways
    a = np.arange(0,1000, dtype='float32')
    a1 = a[np.newaxis, :]
    a1 = a1.repeat(1000, axis=0)
    
    a2 = np.broadcast_to(a, (1000,1000))
    
    # examine size of objects
    sys.getsizeof(a1) # 4000112
    sys.getsizeof(a2) # 112
    
    # try writing
    band.WriteArray(a1) # writes fine
    band.WriteArray(a2) # crash