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))