Different results reading categorical geotiff across platforms

[mjwoods]

When reading land use categories in geotiff format, I get different results on Windows and MacOS binary distributions of R. In both cases, I have tested R_4.5.2 and stars_0.6-8, downloaded as binary packages using Windows 11 and MacOS Tahoe.

The dataset being used was downloaded and unpacked from ABARES catchment scale land use of Australia.

A map showing the expected contents of the dataset is available here.

On each platform, I ran the following code:

library(stars)
library(ggplot2)

tiff <- "clum_50m_2023_v2/clum_50m_2023_v2.tif"

raw <- read_stars(tiff)
rasterio <- list(nXOff=1, nYOff=1,
                 nXSize=dim(raw)[1], nXSize=dim(raw)[2],
                 nBufXSize=dim(raw)[1]/50, nBufYSize=dim(raw)[1]/50,
                 bands=1, resample="nearest_neighbour")

land <- read_stars(tiff, RasterIO=rasterio, RAT="SIMP")
print(ggplot() + geom_stars(data=land))

MacOS results are shown in the first figure below. It appears that numeric values were read instead of categories.

Windows results are as follows. Although categories were read, many categories and data points are missing.

Does this indicate bugs and/or missing features in the binary packages on each platform?

[edzer]

Could you pls give the output of sf::sf_extSoftVersion() and packageVersion("sf") for each platform? sf is where the compiled code is, and the GDAL code linked to.

[mjwoods]

Hi! Thanks for the quick response.

The packages for MacOS are:

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
      "3.13.0"        "3.8.5"        "9.5.1"         "true"         "true"        "9.5.1" 

> packageVersion("sf")
[1] ‘1.0.22’

And for Windows:

> sf::sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H           PROJ 
      "3.13.1"       "3.11.4"        "9.7.0"         "true"         "true"        "9.7.0" 

> packageVersion("sf")
[1] ‘1.0.22’

[kadyb]

Support for *.vat.dbf from ArcGIS was added in GDAL 3.11:

GTiff driver:
add support to read RAT stored in ArcGIS style .tif.vat.dbf auxiliary file

[mjwoods]

I guess that would explain the lack of factor levels on MacOS, where GDAL is an older version. Do you have any idea what is causing the missing values on Windows?

[edzer]

The problem probably is that category labels are repeated. stars molds this into a factor, but R doesn't like duplicate factor levels:

> factor(1:4, levels = c("a", "a", "b", "b"))
Error in `levels<-`(`*tmp*`, value = as.character(levels)) : 
  factor level [2] is duplicated
> structure(1:4, levels = c("a", "a", "b", "b"), class = "factor")
[1] a a b b
Levels: a a b b
Warning message:
In print.factor(x) : duplicated level [2] in factor

if you print the matrix land[[1]] you also get:

Warning message:
In print.factor(x) : duplicated level [2] in factor

and because

> levels(land[[1]]) |> length()
[1] 189
> levels(land[[1]]) |> unique() |> length()
[1] 19

we'd need to merge duplicated levels before constructing something that is will behave like a factor.

[mjwoods]

Thanks for looking into this issue. 😄

I agree that duplicate factor levels would explain problems in the example I showed, but there are other columns in the raster attribute table that do not have duplicate factor levels, and a similar problem occurs in the map (on Windows). For example:

land <- read_stars(tiff, RasterIO=rasterio, RAT="TERTV8")
> levels(land[[1]]) |> unique() |> length()
[1] 189
print(ggplot() + geom_stars(data=land))

The resulting image is shown below, but it has many more missing values than the reference map image, which is shown further down. Does this indicate a bug in the Windows package, or should I be reading the data differently?