Geographic Information Systems Asked on April 8, 2021
I am interested in the terra
R package as it appears to be a successor to the raster
package that is being more actively developed than the raster
package. But it seems like some features are missing.
For example I have an aggregated National Landcover Dataset raster in .GRD format.
If I call x <- raster::raster('nlcd_agg.grd')
the resulting object has a slot x@data@attributes
with the raster attribute table that includes the names of the land cover types corresponding to the integer values in the raster.
However if I call x <- terra::rast('nlcd_agg.grd')
I cannot find that attribute table in the object nor do I see how to associate the object x
with the attribute table of the raster.
I would like to be able to load the raster and have the associated attribute table with land cover type names and default plotting colors load as well. I understand that .GRD is a native format for the raster
package but is there some way to achieve this with the terra
or stars
packages?
Here are two small files that can replicate the behavior (nlcd_agg.grd
and the associated nlcd_agg.gri
):
See the cats
and rats
methods. Or levels
which is the same as cats
.
cats
is for categories (ID and label), rats
for Raster Attribute Tables (ID and n labels). cats are more limited than rats, but generally easier to deal with. They are used to have rasters that behave like "factors".
With the specific file, the problem is that terra
uses GDAL to read it, and the GDAL driver for GRD/GRI is rather incomplete. RATS are read fine, I think, when using e.g. GeoTIFF. I may write a translator to make it easier to extract this from GRD/GRI. But for now, what you can do is something along these lines:
library(raster)
f <- "nlcd_agg.grd"
r <- raster(f)
library(terra)
x <- rast(r)
# get the attributes
lev <- levels(r)[[1]]
Now set the attributes
# this fails because of a bug
rats(x) <- lev
# work around
rat <- lev
rats(x) <- lev
str(rats(x))
They are there, but they are ignored. I would use cats/levels.
lev <- lev[, c("ID", "Land.Cover.Class")]
lev[,2] <- as.character(lev[,2])
x <- rast(r)
levels(x) <- lev
is.factor(x)
#[1] TRUE
x
#class : SpatRaster
#dimensions : 101, 121, 1 (nrow, ncol, nlyr)
#resolution : 3750, 3750 (x, y)
#extent : 1394535, 1848285, 1722765, 2101515 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#source : memory
#name : nlcd_2011_landcover_2011_edition_2014_03_31
#min value : Unclassified
#max value :
# legend shows class names
plot(x)
Labels are returned as cell values
x[c(7357, 5047, 7360, 9307)]
# nlcd_2011_landcover_2011_edition_2014_03_31
#1 Perennial Snow/Ice
#2 Open Water
#3 Perennial Snow/Ice
#4 Unclassified
And the categories are stored if you save the raster to a GeoTIFF
z <- writeRaster(x, "test.tif", overwrite=TRUE)
z
You can also set the color-table
lev <- levels(r)[[1]]
rgb <- lev[,c("Red", "Green", "Blue", "Opacity")]
coltab(x) <- round(rgb*255)
plot(x)
But I did not consider the case where there is a color table and attributes. To be done.
Answered by Robert Hijmans on April 8, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP