rr-sfterrarayshader
Map colours not showing correctly, all SpatRaster values converted to 1 when using crop and mask
I am following this tutorial to create a 3d land cover map and this is the complete code:
library(terra)
library(giscoR)
library(sf)
library(tidyverse)
library(ggtern)
library(elevatr)
library(png)
library(rayshader)
library(magick)
# 2. COUNTRY BORDERS ----
country_sf <- giscoR::gisco_get_countries(
country = "BA",
resolution = "1"
)
plot(sf::st_geometry(country_sf))
png("bih-borders.png")
plot(sf::st_geometry(country_sf))
dev.off()
# 3. DOWNLOAD ESRI LAND COVER TILES ----
urls <- c(
"https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/33T_20220101-20230101.tif",
"https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/34T_20220101-20230101.tif"
)
for(url in urls){
download.file(
url = url,
destfile = basename(url),
mode = "wb"
)
}
# 4. LOAD TILES ----
raster_files <- list.files(
path = getwd(),
pattern = "20230101.tif$",
full.names = T
)
crs <- "EPSG:4326"
for(raster in raster_files){
rasters <- terra::rast(raster)
country <- country_sf |>
sf::st_transform(
crs = terra::crs(
rasters
)
)
land_cover <- terra::crop(
rasters,
terra::vect(
country
),
snap = "in",
mask = T
) |>
terra::aggregate(
fact = 5,
fun = "modal"
) |>
terra::project(crs)
terra::writeRaster(
land_cover,
overwrite=TRUE,
paste0(
raster,
"_bosnia",
".tif"
)
)
}
# 5. LOAD VIRTUAL LAYER ----
r_list <- list.files(
path = getwd(),
pattern = "_bosnia",
full.names = T
)
land_cover_vrt <- terra::vrt(
r_list,
"bosnia_land_cover_vrt.vrt",
overwrite = T
)
# 6. FETCH ORIGINAL COLOURS ----
ras <- terra::rast(
raster_files[[1]]
)
raster_color_table <- do.call(
data.frame,
terra::coltab(ras)
)
head(raster_color_table)
hex_code <- ggtern::rgb2hex(
r = raster_color_table[,2],
g = raster_color_table[,3],
b = raster_color_table[,4]
)
# 7. ASSIGN COLORS TO RASTER----
cols <- hex_code[c(2:3, 5:6, 8:12)]
from <- unique_vals
to <- t(col2rgb(cols))
land_cover_vrt <- na.omit(land_cover_vrt)
land_cover_bosnia <- terra::subst(
land_cover_vrt,
from = from,
to = to,
names = cols
)
terra::plotRGB(land_cover_bosnia)
# 8. DIGITAL ELEVATION MODEL ----
elev <- elevatr::get_elev_raster(
locations = country_sf,
z = 9, clip = "locations"
)
crs_lambert <-
"+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_frfs"
land_cover_bosnia_resampled <- terra::resample(
x = land_cover_bosnia,
y = terra::rast(elev),
method = "near"
) |>
terra::project(crs_lambert)
terra::plotRGB(land_cover_bosnia_resampled)
img_file <- "land_cover_bosnia.png"
terra::writeRaster(
land_cover_bosnia_resampled,
img_file,
overwrite = T,
NAflag = 255
)
img <- png::readPNG(img_file)
# 9. RENDER SCENE ----
elev_lambert <- elev |>
terra::rast() |>
terra::project(crs_lambert)
elmat <- rayshader::raster_to_matrix(
elev_lambert
)
h <- nrow(elev_lambert)
w <- ncol(elev_lambert)
elmat |>
rayshader::height_shade(
texture = colorRampPalette(
cols[9]
)(256)
) |>
rayshader::add_overlay(
img,
alphalayer = 1
) |>
rayshader::plot_3d(
elmat,
zscale = 12,
solid = F,
shadow = T,
shadow_darkness = 1,
background = "white",
windowsize = c(
w / 5, h / 5
),
zoom = .5,
phi = 85,
theta = 0
)
rayshader::render_camera(
zoom = .58
)
# 10. RENDER OBJECT ----
u <- "https://dl.polyhaven.org/file/ph-assets/HDRIs/hdr/4k/air_museum_playground_4k.hdr"
hdri_file <- basename(u)
download.file(
url = u,
destfile = hdri_file,
mode = "wb"
)
filename <- "3d_land_cover_bosnia-dark.png"
rayshader::render_highquality(
filename = filename,
preview = F,
light = F,
environment_light = hdri_file,
intensity_env = 1,
rotate_env = 90,
interactive = F,
parallel = T,
width = w * 1.5,
height = h * 1.5
)
# 11. PUT EVERYTHING TOGETHER ----
c(
"#419bdf", "#397d49", "#7a87c6",
"#e49635", "#c4281b", "#a59b8f",
"#a8ebff", "#616161", "#e3e2c3"
)
legend_name <- "land_cover_legend.png"
png(legend_name)
par(family = "mono")
plot(
NULL,
xaxt = "n",
yaxt = "n",
bty = "n",
ylab = "",
xlab = "",
xlim = 0:1,
ylim = 0:1,
xaxs = "i",
yaxs = "i"
)
legend(
"center",
legend = c(
"Water",
"Trees",
"Crops",
"Built area",
"Rangeland"
),
pch = 15,
cex = 2,
pt.cex = 1,
bty = "n",
col = c(cols[c(1:2, 4:5, 9)]),
fill = c(cols[c(1:2, 4:5, 9)]),
border = "grey20"
)
dev.off()
# filename <- "land-cover-bih-3d-b.png"
lc_img <- magick::image_read(
filename
)
my_legend <- magick::image_read(
legend_name
)
my_legend_scaled <- magick::image_scale(
magick::image_background(
my_legend, "none"
), 2500
)
p <- magick::image_composite(
magick::image_scale(
lc_img, "x7000"
),
my_legend_scaled,
gravity = "southwest",
offset = "+100+0"
)
magick::image_write(
p, "3d_bosnia_land_cover_final.png"
)
I am using MacOS, R version = 4.4.0, rayshader version = 0.35.7, terra version = 1.7-78.
This is my output:
However, this is what the output should be like:
I am using MacOS, R version = 4.4.0, rayshader version = 0.35.7, terra version = 1.7-78.
Solution
For some reason, the values in large SpatRaster files are all converted to 1 if you call mask inside crop(). The solution is to run mask() separately. I'm not sure whether this is caused by a hardware issue or not, but if you replace the for() loop in step 4 with:
for(raster in raster_files){
rasters <- rast(raster)
country <- country_sf |>
st_transform(crs(rasters))
land_cover <- crop(rasters, vect(country), snap = "in")
land_cover <- mask(land_cover, vect(country)) |>
aggregate(fact = 5, fun = "modal") |>
project(crs)
writeRaster(land_cover,
paste0(raster, "_bosnia", ".tif"),
overwrite = TRUE)
}
then run the rest of the code, the output 3d_bosnia_land_cover_final.png file will show the correct colours.