array orientation - quick explainer

[mdsumner]

We recommend using "raster order", that is how GDAL and packages terra, stars, gdalraster, and vapour deal with pixel values - this is from the topleft of the map, to the right across the top row and then each row after that. A matrix oriented this way is what rasterImage() expects, which is different from what image() expects. Of course we have to deal with the existing convention in R which is about how matrix() populates and the way we get values from a matrix if we flatten it.

Note that the spatial packages in R all have this covered in their own plotting mechs, but if you want to be able to navigate the underlying details, read on. (Note that we aren't considering "tiles" here, if your source data is tiled the native storage is this within each separate tile (aka "chunk"), so this still applies at that lower level of organization. With GDAL we can read a native tile or across them however we like, we're ignoring that detail here).

So 3 things

• raster order is topleft to right and down so for matrix(, nrow, byrow = TRUE) is how to generate a matrix in raster order from GDAL values
• for image() to use that it must flip and transpose t(x[nrow(x):1, ])
• if we want the values in the original order from a matrix created in that orientation we must t() it

 ## ximage has an example matrix in "raster"/rasterImage orientation
 x <- ximage::topo
 
 ## or could be this as equivalent
 #r <- terra::project(terra::rast("/vsicurl/https://projects.pawsey.org.au/idea-gebco-tif/GEBCO_2024.tif"), 
 #                                terra::rast(terra::ext(-180, 180, -90, 90), ncols = 720, nrows = 360, crs = "EPSG:4326"), 
 #                  by_util = TRUE)

 #x <- matrix(terra::values(r)[,1], nrow(r), byrow = TRUE)


 dim(t(x[nrow(x):1, ]))
#> [1] 720 360

 ## Antarctica is at the end, we transpose else we go down the rows first rather then left-right
 plot(c(t(x)), pch = ".")

 image(x)  ## nope

 image(t(x[nrow(x):1, ]))  ## yep

 plot(c(-180, 180), c(-90, 90))
 rasterImage(scales::rescale(x), -180, -90, 180, 90)  ## yep

 ## so ximage

 ximage::ximage(x)  ## index space

 ximage::ximage(x, c(-180, 180, -90, 90))  ## geog space

^{Created on 2025-04-11 with reprex v2.1.1}