Search code examples
rggplot2rasterartifactselevatr

Gridded dot artifacts in geom_raster plot


I'm attempting to make a high resolution figure of a digital elevation model (raster), but my plots appear to have some distorted features. Here's a reproducible example:

library(elevatr)
library(raster)
library(tidyverse)
library(sp)

ext <- extent(-864434.2, -771071.7, 4019064, 4054974) %>%
  as("SpatialPolygons") %>%
  spsample(100, "regular")

crs(ext) <- "+proj=utm +zone=19 +ellps=GRS80 +datum=NAD83"

alt <- get_elev_raster(ext, 
                       prj = crs(ext), 
                       z = 11) %>%
  crop(ext)
names(alt) <- "alt"

ggplot() + 
  geom_raster(data = as.data.frame(alt, xy = T), aes(x = x, y = y, fill = alt)) +
  scale_fill_gradientn("Elevation (m)", colors = terrain.colors(256), na.value = NA) +
  theme(legend.position = "none")
ggsave("gg_gsmnp_alt_test.jpeg", width = 8, height = 4, units = "in", dpi = 600)

Plot with artifact #1 Plot with artifact #2

The artifacts become more prevalent with increased resolution (if I set the z parameter to a higher resolution, e.g. 11 or 12, the artifacts are much more abundant).

Calculating a hill shade layer exacerbates these artifacts:

hs <- hillShade(terrain(alt, opt = "slope"),
                terrain(alt, opt = "aspect"))

ggplot() + 
  geom_raster(data = as.data.frame(alt, xy = T), aes(x = x, y = y, fill = alt)) +
  geom_raster(data = as.data.frame(hs, xy = T), aes(x = x, y = y, alpha = 1 - layer), fill = "gray20") +
  scale_fill_gradientn("Elevation (m)", colors = terrain.colors(256), na.value = NA) +
  scale_alpha(guide = FALSE, range = c(0,1))
ggsave("gg_gsmnp_alt_hs_test.jpeg", width = 8, height = 4, units = "in", dpi = 600)

Plot with artifact #3

Furthermore, if I aggregate the layer using bilinear interpolation, the artifacts predictably disappear:

alt_agg <- alt %>%
   aggregate(fact = 2)
hs_agg <- hs %>%
  aggregate(fact = 2)

ggplot() + 
  geom_raster(data = as.data.frame(alt_agg, xy = T), aes(x = x, y = y, fill = alt)) +
  geom_raster(data = as.data.frame(hs_agg, xy = T), aes(x = x, y = y, alpha = 1 - layer), fill = "gray20") +
  scale_fill_gradientn("Elevation (m)", colors = terrain.colors(256), na.value = NA) +
  scale_alpha(guide = FALSE, range = c(0,1)) +
  theme(legend.position = "none")
ggsave("gg_gsmnp_alt_hs_agg_test.jpeg", width = 8, height = 4, units = "in", dpi = 600)

Plot with interpolated values (no artifacts)

While the graphical issues are minor, I'm curious if these artifacts may contribute to anomalies in data extracted from the raster. For instance, if I extract data using raster::extract at these locations, will I receive spurious values?

Any advice or answers would be greatly appreciated.

Best,

-Alex. (he/him/his)


Solution

  • Better late than never... Apologies.

    Root cause of this was GDAL warp when I tried to transform and moscaic at the same time. I now use two steps, mosaic, then transform. Version at https://github.com/jhollist/elevatr has the fix. I am working on a CRAN release. Hope this helps and again apologies.