White dots on SpatRaster when plotting with tidyterra::geom_spatraster_RGB() and not terra::plotRGB()?

I'm making a map of this island and whenever I try to plot it using tidyterra::geom_spatraster_RGB() there are these white dots and the graph is more pale than when I use terra::plot_RGB(). I definitely want to use ggplot to plot this data because I have points that I want to use that are in another dataframe and I use geom_point to plot them. I also eventually want to apply ggmagnify (geom_magnify()) to make a map inset to zoom into a specific area on the island (which I have gotten to work but the only issue is the little white dots!!).

My goal output is a SpatRaster or something I can use in ggplot2 that doesn't have the white dots or clouds obstructing the geology/ coastline of the island. (Basically the terra::plot_RGB() plot except useable in ggplot2).

I'm pretty new to GIS in general and especially GIS in R, thanks in advance for any help!

library(rsi)
library(tidyterra)
library(terra)
library(tidyverse)
library(sf)

annette_boundary = st_point(c(-131.5, 55.15))
annette_boundary = st_sfc(annette_boundary, crs=4326)
annette_boundary = st_buffer(st_transform(annette_boundary, 4326), 20000)

ab_sentinel2 = get_sentinel2_imagery(
  annette_boundary,
  start_date = "2024-10-20",
  end_date = "2024-10-22",
  output_filename = "sentinel2_imagery.tif"
  
)


ab_sentinel2_rast = terra::rast("sentinel2_imagery.tif")
ab_sentinel2_rast = terra::project(ab_sentinel2_rast, "EPSG:4326")

# plot with ggplot
ggplot(data = ab_sentinel2_rast)+ 
  tidyterra::geom_spatraster_rgb(data = ab_sentinel2_rast, r=4, g=3, b=2, 
                                 interpolate=F, max_col_value = 5000, maxcell = 5e+05)

# plot with terra
terra::plotRGB(ab_sentinel2_rast, r=4, g=3, b=2, stretch="lin")

The white dots are caused by lack of pixels (colors/values) in a original image.

First of all, take the image (please note I had to convert the 4326 to projected coordinates (EPSG:6394), otherwise wasn't able to download the image using {rsi} function:

library(rsi)
library(tidyterra)
library(terra)
library(tidyverse)
library(sf)
sf::sf_use_s2(TRUE)
annette_boundary <- st_point(c(-131.5, 55.15))
annette_boundary <- st_sfc(annette_boundary, crs="EPSG:4326") |>
  sf::st_as_sf()

annette_boundary <- annette_boundary |>
  sf::st_transform(crs = "EPSG:6394") |>
  st_buffer(dist = 20000) |>
  sf::st_sf()

ab_sentinel2 = get_sentinel2_imagery(
  annette_boundary,
  pixel_x_size = 10,
  pixel_y_size = 10,
  start_date = "2024-10-20",
  end_date = "2024-10-22",
  output_filename = "sentinel2_imagery.tif"
)

ab_sentinel2_rast <- terra::rast("sentinel2_imagery.tif")
ab_sentinel2_rast

#> class       : SpatRaster 
#> dimensions  : 4000, 4000, 12  (nrow, ncol, nlyr)
#> resolution  : 10, 10  (x, y)
#> extent      : 936802.9, 976802.9, 351310.2, 391310.2  (xmin, xmax, ymin, ymax)
#> coord. ref. : NAD83(2011) / Alaska zone 1 (EPSG:6394) 
#> source      : sentinel2_imagery.tif 
#> names       : A, B, G, R, RE1, RE2, ... 

Now, let's plot it:

terra::plotRGB(ab_sentinel2_rast, r=4, g=3, b=2, stretch="lin")

The white dots are visible, but let's check those layers:

r <- terra::values(ab_sentinel2_rast[[4]])
g <- terra::values(ab_sentinel2_rast[[3]])
b <- terra::values(ab_sentinel2_rast[[2]])
which(is.na(r) & is.na(g) & is.na(b))
#>   [1]    143    144    209    224    225    226    227    307   1480   1608   2403   2416   2418   2426   2427
#>  [16]   2919   2927   3003   3021   3022   3025   3223   3970   3971   4142   4151   4152   4153   4276   4277
[...]

Let's "repair" our raster substituting the missing values with mean from neigbourghs:

ab_sentinel2_rast[[4]] <- terra::focal(ab_sentinel2_rast[[4]], fun = "mean", na.policy = "only")
ab_sentinel2_rast[[3]] <- terra::focal(ab_sentinel2_rast[[3]], fun = "mean", na.policy = "only")
ab_sentinel2_rast[[2]] <- terra::focal(ab_sentinel2_rast[[2]], fun = "mean", na.policy = "only")

and plot it:

ggplot(data = ab_sentinel2_rast)+ 
  tidyterra::geom_spatraster_rgb(data = ab_sentinel2_rast, r=4, g=3, b=2, 
                                 interpolate=F, max_col_value = 5000, maxcell = 5e+05)

There is probably still few dots, you can increase the number of window in focal() function with w=5 or similar.

You can add your points or other features to plot made with terra::plot() as well, have a look on points(), lines() and polys() function in {terra} package.