1

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: this

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

0

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
  • Thank you so much!!! I followed your sample code and managed to make it (with 1 problem...) I did have to change: xmin,ymin,xmax,ymax = [lons.min()-xres/2, lats.min()-yres/2,lons.max()+xres/2, lats.max()+yres/2] since the coords were centroids; and array[y, x] = preds[j, i] since the shape was (94,1461). – plummms Apr 8 at 6:01
  • The 1 problem: The final raster is inverted on the y-axis – plummms Apr 8 at 6:02
  • Solved my problem by adding np.flipud(array) before the WriteArray(array) line. Source – plummms Apr 8 at 6:31

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.