Sunday, 21 August 2016

How to write a single pixel or a line of pixels to GeoTIFF using GDAL/Python


Let's say I want to write a single pixel or a 1D array of pixels (e.g.all at the same latitude) of pixels to a single-band of a Float32 GeoTIFF using the Python GDAL bindings. How can I do so? Or is it not possible?


I currently have the following method for truly 2D ndarrays. But occasionally I get requests for a pixel or line of pixels, which cause this method to fail.



raster_band.WriteArray(my_pixel_array)

FURTHER CLARIFICATION:


I'm trying to write a new image where the entire extent of the new image is either just one pixel or just a line of pixels. Thus, I've tried just setting xoff and yoff to 0, and this fails the same way as without xoff and yoff. The error I'm getting is expected array of dim 2, e.g. when the shape of my ndarray is (4,) (and when the GDAL raster is initialized with RasterXSize to 4 and RasterYSize to 1 via a pertinent out_raster.SetGeoTransform (EPSG:4326 location and size) method).



Answer



The WriteArray() method has two useful optional arguments to accomplish this:



  • xoff: specifies the offset on the x axis.

  • yoff: specifies the offset on the y axis.



This means that it will write a 2D array (this is very important) to the band, starting from the specified offsets (which are 0 by default). This is very useful because you don't have to read the whole array into memory, change one pixel or one row, and then write it again. This lets you specify where you want to write the data in the array. Imagine you wanted to replace all the pixels in the last row of the raster with values of 100, and you also wanted to replace the most upper right pixel with a value of 99.


import gdal
import numpy as np

# open the raster in writing mode
ds = gdal.Open('path/to/raster.tif', 1)

# get first band and raster dimensions
band = ds.GetRasterBand(1)
cols = ds.RasterXSize

rows = ds.RasterYSize

# create the arrays you are going to write
row_data = np.full(shape=(1, cols), fill_value=100) # note that eventhough it contains just one row of data it is a 2D array
pixel_data = np.array([[99]]) # again, a 2D array. This time it contains just one value

# write the arrays specifying the proper offsets
band.WriteArray(row_data, xoff=0, yoff=(rows-1))
band.WriteArray(pixel_data, xoff=(cols-1), yoff=0)


# save and close the file
del band, ds

Note that the WriteArray() method only accepts a 2D array, whether it has just one value, one row, or many rows. Furthermore, if you try to write an array with higher dimensions (e.g. you set the offset to be the last row but pass an array with two rows or more columns that ds.RasterXSize), it will throw an Exception because the array does not fit.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...