Writing numpy arrays to irregularly shaped, multiband raster

Question

I have the following numpy arrays:

predictions; type is float32; shape: (94, 1461)

lats; float32; shape: (94,); centroid of cell in vector layer, decimal degs.

lons; float32; shape: (94,); centroid of cell in vector layer, decimal degs.

• 94 corresponds to the number of cells in the vector grid below
• 1461 corresponds to the number of bands

I've looked at similar questions but most of them discuss exporting full (i.e. no missing, rectangular) rasters or single-band rasters. I was wondering how could I export a multi-band raster with multiple noData cells similar to the gridded vector layer here:

Additional info: spatial ref: epsg:4326, GTiff, cell size 0.0270

My attempt:

xmin,ymin,xmax,ymax = [lons.min(),lats.min(),lons.max(),lats.max()]
xres = 0.0270
yres = 0.0270
geotransform=(xmin,xres,0,ymax,0, -yres)
output_raster = gdal.GetDriverByName('GTiff').Create(output_file, 9, 17, 1461, gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)  
for i in range(1461):
        output_raster.GetRasterBand(i+1).WriteArray(predictions)
        output_raster.GetRasterBand(i+1).SetNoDataValue(noData)
output_raster.FlushCache()

Throws a: line 229, output_raster.GetRasterBand(i+1).WriteArray(predictions) ValueError: array larger than output file, or offset off edge

Answer 1

You can't, you have to generate a regularly shaped array.

i.e.

import math

from osgeo import gdal
import numpy as np

xres = 0.0270
yres = 0.0270
nrows = 17
ncols = 9
nbands = 1461
nodata = -9999
ncells = 94

# ======= Make some dummy data =======
output_file = '/tmp/test.tif'
xmin, ymin = 120, -30
cells = np.random.choice(np.arange(nrows*ncols), ncells, replace=False)

lats = np.arange(ymin, ymin+nrows*yres, yres)
lons = np.arange(xmin, xmin+ncols*xres, xres)
lats, lons = np.meshgrid(lats, lons)
lats = lats.ravel()[cells]
lons = lons.ravel()[cells]

predictions = np.random.random((nbands, ncells))
# ======= End make some dummy data =======

# Make an empty 1 band array to fill with predictions
array = np.empty((nrows, ncols), dtype=np.float32)

xmin,ymin,xmax,ymax = [lons.min(), lats.min(),lons.max(), lats.max()]
geotransform=(xmin, xres, 0, ymax, 0, -yres)

output_raster = gdal.GetDriverByName('GTiff').Create(output_file, ncols, nrows, nbands, gdal.GDT_Float32)  # Open the file
output_raster.SetGeoTransform(geotransform)

# Loop bands
for i in range(nbands):

    # Init array with nodata
    array[:] = nodata

    # Loop lat/lons inc. index j
    for j, (lon, lat) in enumerate(zip(lons, lats)):
        # Calc x, y pixel index
        x = math.floor((lon - xmin) / xres)
        y = math.floor((lat - ymin) / xres)

        # Fill the array at y, x with the value from predictions at band i, index j 
        array[y, x] = predictions[i, j]

    output_raster.GetRasterBand(i+1).WriteArray(array)
    output_raster.GetRasterBand(i+1).SetNoDataValue(nodata)

del output_raster