rggplot2terratidyterra
issues plotting cropped raster DEM in ggplot
It's probably a stupid mistake but I can't figure out what I am doing wrong. I have a DEM raster I can visualize using plot() and ggplot() just fine. But when I try to crop and mask it to a smaller polygon, it introduces NAs such that calling the cropped raster object shows the min/max values as NA.
This causes a problem when plotting using ggplot()+ geom_spatraster() as it gives just three values to be represented and NA, and not the continuous range of values in the in the actual raster. my code and data are below:
#The mask object, called dogwood.sf in the code
structure(list(Array = "Dogwood", bbox = structure(list(structure(list(
structure(c(456582.484888111, 459612.80088017, 459612.80088017,
456582.484888111, 456582.484888111, 4038674.29664676, 4038674.29664676,
4041198.55641265, 4041198.55641265, 4038674.29664676), dim = c(5L,
2L))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON",
"sfc"), precision = 0, bbox = structure(c(xmin = 456582.484888111,
ymin = 4038674.29664676, xmax = 459612.80088017, ymax = 4041198.55641265
), class = "bbox"), crs = structure(list(input = "EPSG:26915",
wkt = "PROJCRS[\"NAD83 / UTM zone 15N\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"UTM zone 15N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-93,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"North America - between 96°W and 90°W - onshore and offshore. Canada - Manitoba; Nunavut; Ontario. United States (USA) - Arkansas; Illinois; Iowa; Kansas; Louisiana; Michigan; Minnesota; Mississippi; Missouri; Nebraska; Oklahoma; Tennessee; Texas; Wisconsin.\"],\n BBOX[25.61,-96,84,-90]],\n ID[\"EPSG\",26915]]"), class = "crs"), n_empty = 0L)), row.names = c(NA,
-1L), sf_column = "bbox", agr = structure(c(Array = NA_integer_), levels = c("constant",
"aggregate", "identity"), class = "factor"), class = c("sf",
"tbl_df", "tbl", "data.frame"))
#not too sure how to share the DEM, but in this example, it's called DEM
terra::crop(DEM,dogwood.sf)->dem.crop
terra::mask(dem.crop,dogwood.sf)->dem.mask
terra::plot(dem.mask) #works as expected
ggplot2::ggplot()+
tidyterra::geom_spatraster(dem.mask) # does not work, showing only 3 values
dem.mask # see the min/max values are NA
terra::values(dem.mask) # see full continuous range of raster values
ggplot2::ggplot()+
tidyterra::geom_spatraster(dem.mask)+
tidyterra::scale_fill_terrain_c()
# gives error about supplying discrete values to a continuous scale.
I appreciate any help anyone cares to provide.
Solution
I could not replicate your issue using a sample DEM, so there's potentially something going on in your session environment or some other workflow issue.
Whatever the cause, here's a reprex to see if your issue persists using a different DEM source. At the very least, this may help you to debug your issue. Note that both datasets default to the same CRS e.g. NAD83 UTM Zone 15/EPSG:26915.
library(sf)
library(terra)
library(elevatr) # For downloading DEM
library(ggplot2)
library(tidyterra)
# Load sf
dogwood.sf <- structure(list(Array = "Dogwood", bbox = structure(list(structure(list(
structure(c(456582.484888111, 459612.80088017, 459612.80088017,
456582.484888111, 456582.484888111, 4038674.29664676, 4038674.29664676,
4041198.55641265, 4041198.55641265, 4038674.29664676), dim = c(5L,
2L))), class = c("XY", "POLYGON", "sfg"))), class = c("sfc_POLYGON",
"sfc"), precision = 0, bbox = structure(c(xmin = 456582.484888111,
ymin = 4038674.29664676, xmax = 459612.80088017, ymax = 4041198.55641265
), class = "bbox"), crs = structure(list(input = "EPSG:26915",
wkt = "PROJCRS[\"NAD83 / UTM zone 15N\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"UTM zone 15N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-93,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"North America - between 96°W and 90°W - onshore and offshore. Canada - Manitoba; Nunavut; Ontario. United States (USA) - Arkansas; Illinois; Iowa; Kansas; Louisiana; Michigan; Minnesota; Mississippi; Missouri; Nebraska; Oklahoma; Tennessee; Texas; Wisconsin.\"],\n BBOX[25.61,-96,84,-90]],\n ID[\"EPSG\",26915]]"), class = "crs"), n_empty = 0L)), row.names = c(NA,
-1L), sf_column = "bbox", agr = structure(c(Array = NA_integer_), levels = c("constant",
"aggregate", "identity"), class = "factor"), class = c("sf",
"tbl_df", "tbl", "data.frame"))
# Download example raster DEM tile
DEM <- get_elev_raster(dogwood.sf, z = 14, expand = 1)
# Convert DEM to SpatRaster and assign values
DEM <- rast(DEM)
DEM[] <- values(DEM)
# Crop and mask DEM
dem.mask <- crop(DEM, dogwood.sf, mask = TRUE)
dem.mask
# class : SpatRaster
# dimensions : 658, 790, 1 (nrow, ncol, nlyr)
# resolution : 3.834583, 3.834583 (x, y)
# extent : 456583.7, 459613, 4038676, 4041199 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs
# source(s) : memory
# varname : fileb6468fc3331
# name : fileb6468fc3331
# min value : 302.2790
# max value : 430.6848
ggplot() +
geom_spatraster(data = dem.mask) +
geom_sf(data = dogwood.sf, fill = NA, colour = "#DF536B", linewidth = .5) +
coord_sf(datum = 26915) +
scale_x_continuous(breaks = seq(from = st_bbox(dogwood.sf)[[1]],
to = st_bbox(dogwood.sf)[[3]], length.out = 4))
