Raster Data

In this section, we revisit the gdalworkshop/world.tif dataset.

dataset = AG.read("gdalworkshop/world.tif")

GDAL Dataset (Driver: GTiff/GeoTIFF)
File(s): 

Dataset (width x height): 2048 x 1024 (pixels)
Number of raster bands: 3
  [GA_ReadOnly] Band 1 (Red): 2048 x 1024 (UInt8)
  [GA_ReadOnly] Band 2 (Green): 2048 x 1024 (UInt8)
  [GA_ReadOnly] Band 3 (Blue): 2048 x 1024 (UInt8)

A description of the display is available in Raster Datasets.

Raster Bands

We can examine an individual raster band

band = ArchGDAL.getband(dataset, 1)

[GA_ReadOnly] Band 1 (Red): 2048 x 1024 (UInt8)
    blocksize: 256×256, nodata: -1.0e10, units: 1.0px + 0.0
    overviews: (0) 1024x512 (1) 512x256 (2) 256x128 
               (3) 128x64 (4) 64x32 (5) 32x16 
               (6) 16x8 

You can programmatically retrieve the information in the header using

• ArchGDAL.accessflag(band): the access flag for this band. (GA_ReadOnly)
• colorinterp =ArchGDAL.getcolorinterp(band): color interpretation of the values in the band (GCI_RedBand)
• ArchGDAL.getname(colorinterp): name (string) corresponding to color interpretation ("Red")
• ArchGDAL.width(band): width (pixels) of the band (2048)
• ArchGDAL.height(band): height (pixels) of the band (1024)
• ArchGDAL.indexof(band): the index of the band (1+) within its dataset, or 0 if unknown. (1)
• ArchGDAL.pixeltype(band): pixel data type for this band. (UInt8)

You can get additional attribute information using

• ArchGDAL.getscale(band): the scale in units = (px * scale) + offset (1.0)
• ArchGDAL.getoffset(band): the offset in units = (px * scale) + offset (0.0)
• ArchGDAL.getunittype(band): name for the units, e.g. "m" (meters) or "ft" (feet). ("")
• ArchGDAL.getnodatavalue(band): a special marker value used to mark pixels that are not valid data. (-1.0e10)
• (x,y) =ArchGDAL.blocksize(band): the "natural" block size of this band ((256,256))

Finally, you can obtain overviews:

• ArchGDAL.noverview(band): the number of overview layers available, zero if none. (7)
• ArchGDAL.getoverview(band, i): returns the i-th overview in the raster band. Each overview is itself a raster band, e.g.

ArchGDAL.getoverview(band, 2)

[GA_ReadOnly] Band 1 (Red): 256 x 128 (UInt8)
    blocksize: 128×128, nodata: -1.0e10, units: 1.0px + 0.0
    overviews: 

Raster I/O

Reading Raster Values

The general operative method for reading in raster values from a dataset or band is to use ArchGDAL.read().

• ArchGDAL.read(dataset): reads the entire dataset as a single multidimensional array.
• ArchGDAL.read(dataset, indices): reads the raster bands at the indices (in that order) into a multidimensional array.
• ArchGDAL.read(dataset, i): reads the i-th raster band into an array.
• ArchGDAL.read(band): reads the raster band into an array.

You can also specify the subset of rows and columns (provided as UnitRanges like 2:9) to read:

• ArchGDAL.read(dataset, indices, rows, cols)
• ArchGDAL.read(dataset, i, rows, cols)
• ArchGDAL.read(band, rows, cols)

On other occasions, it might be easier to first specify a position (xoffset,yoffset) to read from, and the size (xsize, ysize) of the window to read:

• ArchGDAL.read(dataset, indices, xoffset, yoffset, xsize, ysize)
• ArchGDAL.read(dataset, i, xoffset, yoffset, xsize, ysize)
• ArchGDAL.read(band, xoffset, yoffset, xsize, ysize)

You might have an existing buffer that you wish to read the values into. In such cases, the general API for doing so is to write ArchGDAL.read!(source, buffer, args...) instead of ArchGDAL.read(source, args...).

Writing Raster Values

For writing values from a buffer to a raster dataset or band, the following methods are available:

• ArchGDAL.write!(band, buffer)
• ArchGDAL.write!(band, buffer, rows, cols)
• ArchGDAL.write!(band, buffer, xoffset, yoffset, xsize, ysize)
• ArchGDAL.write!(dataset, buffer, i)
• ArchGDAL.write!(dataset, buffer, i, rows, cols)
• ArchGDAL.write!(dataset, buffer, i, xoffset, yoffset, xsize, ysize)
• ArchGDAL.write!(dataset, buffer, indices)
• ArchGDAL.write!(dataset, buffer, indices, rows, cols)
• ArchGDAL.write!(dataset, buffer, indices, xoffset, yoffset, xsize, ysize)

Windowed Reads and Writes

Following the description in mapbox/rasterio's documentation, a window is a view onto a rectangular subset of a raster dataset. This is useful when you want to work on rasters that are larger than your computers RAM or process chunks of large rasters in parallel.

For that purpose, we have a method called ArchGDAL.windows(band) which iterates over the windows of a raster band, returning the indices corresponding to the rasterblocks within that raster band for efficiency:

using Base.Iterators: take  # to prevent showing all blocks
windows = ArchGDAL.windows(band)

for (cols, rows) in take(windows, 5)
    @info "Window" cols rows
end

┌ Info: Window
│   cols = 1:256
└   rows = 1:256
┌ Info: Window
│   cols = 1:256
└   rows = 257:512
┌ Info: Window
│   cols = 1:256
└   rows = 513:768
┌ Info: Window
│   cols = 1:256
└   rows = 769:1024
┌ Info: Window
│   cols = 257:512
└   rows = 1:256

Alternatively, we have another method called ArchGDAL.blocks(band) which iterates over the windows of a raster band, returning the offset and size corresponding to the rasterblocks within that raster band for efficiency:

blocks = ArchGDAL.blocks(band)
for (xyoffset, xysize) in take(blocks, 5)
    @info "Window offset" xyoffset xysize
end

┌ Info: Window offset
│   xyoffset = (0, 0)
└   xysize = (256, 256)
┌ Info: Window offset
│   xyoffset = (1, 0)
└   xysize = (256, 256)
┌ Info: Window offset
│   xyoffset = (2, 0)
└   xysize = (256, 256)
┌ Info: Window offset
│   xyoffset = (3, 0)
└   xysize = (256, 256)
┌ Info: Window offset
│   xyoffset = (0, 1)
└   xysize = (256, 256)