I have been using terra to read and then save various global raster data sets. Today I noticed quite different value ranges depending on how I look at the data or what I do with it. I am confused and nervous to undertake any analyses until I can explain why I am seeing what I am seeing. I would be most grateful if any among you could let me know why the value ranges I am seeing are so different.
Using global road density as an example. I import the GRIP Total density, all types combined ascii raster using terra as follows:
roadDir <- c("downloads/GRIP4_density_total/")
tmp <- paste0(roadDir, "grip4_total_dens_m_km2.asc")
roadDensity <- terra::rast(tmp)
If I then look at roadDensity I see that the range of values is 0 to 99445 (m per square km):
> roadDensity
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source : grip4_total_dens_m_km2.asc
name : grip4_total_dens_m_km2
min value : 0
max value : 99445
If I import the original data into QGIS I get the same cell size, dimensions and value range as shown in the code snippet above so I take those values as the true values in the data. So far so good.
Problems emerge when: a) I look at summaries of the data or b) write the resulting raster to a .tif file for use in further processing.
a) Summary
terra::global(roadDensity, c( "min", "max"), na.rm=TRUE)
min max
grip4_total_dens_m_km2 0 330306
I do not understand why I am now getting a very different max value (330306).
b) Write to .tif
x <- writeRaster(roadDensity,"test.tif", datatype="INT4U", overwrite=TRUE)
x
class : SpatRaster
dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
resolution : 0.08333333, 0.08333333 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : test.tif
name : grip4_total_dens_m_km2
min value : 0
max value : 330306
So terra's global and writeRaster produce the same results.
I understand that through changing the datatype value to one of its possibilities ("INT1U", "INT2U", "INT4U") writeRaster changes the value range because of the max possible values that particular data structure can hold. What I do not understand is how "INT4U" can increase the value range so dramatically beyond the original without any change in cell size. In addition the INT*U value should not change the summary values produced by global.
I am sure its a very simple thing and I am due an "ah duh" moment, but I would appreciate any guidance to aid my understanding and thence analytical confidence.
Thank you
I do not see that. My best guess is that QGIS has created a sidecar file with approximate extreme values.
Get the data in a clean folder
dir.create("test")
setwd("test")
url <- "https://dataportaal.pbl.nl/downloads/GRIP4/GRIP4_density_total.zip"
zip <- basename(url)
if (!file.exists(zip)) {
download.file(url, zip, mode="wb")
unzip(zip)
}
min and max are not reported for an ascii file (as these are not stored in the file format)
library(terra)
#terra 1.7.78
r <- rast("grip4_total_dens_m_km2.asc")
r
#class : SpatRaster
#dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
#resolution : 0.08333333, 0.08333333 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#source : grip4_total_dens_m_km2.asc
#name : grip4_total_dens_m_km2
Make a new SpatRaster to reveal the min and max values:
r * 1
#class : SpatRaster
#dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
#resolution : 0.08333333, 0.08333333 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
#source(s) : memory
#varname : grip4_total_dens_m_km2
#name : grip4_total_dens_m_km2
#min value : 0
#max value : 330306
Or compute them like this:
minmax(r, compute=TRUE)
# grip4_total_dens_m_km2
#min 0
#max 330306
Perhaps you can try this again in an empty folder. I am guessing that in your folder there is a sidecar file left from previous work (e.g. grip4_total_dens_m_km2.xml) that provides these numbers. Perhaps this file was created by QGIS, possibly using a quick estimate. Also please report the version of terra you are using.