Plot Satellite data (.nc format) in ggplot with higher resolution and adding colour bar in R

Question

I would like to plot satellite chlorophyll of a region set by coordinates. However the resolution seems to be very low, (a lot of big blocks/pixels). And

I would like to:

• Extrapolate or smooth the data so the blocks are smaller.
• Plot values in log scale because the gradient is very small (0.01-0.5 units)
• Make the colour bar appear, which is not.

I am using the raster package following this tutorial although they have alternative packages: ncdf4 and ocean map

The nc file (it's 16kb) can be found here

require(oceanmap)
require(raster)
require(ncdf4)
require(tidyverse)
library(sf)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

# Plot of the desired area without chl
  ggplot(data = world) +
  geom_sf() +
  coord_sf(xlim = c(-22, -12), ylim = c(24, 29), expand = FALSE)

 
chl.a = raster::raster("~/Downloads/file_chl_satellite.nc", 
                       varname = "chl")  ## read as raster object

rufiji.extent = raster::extent(-21.75, -13, 24, 29) ## create chopping object

chl.a = raster::crop(x = chl.a, y = rufiji.extent) ## chop the raster object

# assign all cells with values between 0.01 and 0.5
chl.a = chl.a %>% 
  raster::clamp(upper = 0.5, # upper clamp value
                lower = 0.01,  # lower clamp value
                usevalues = FALSE # force value to NA
  )

## convert from raster object to data frame and then tibble
chl.a.tb = chl.a %>% 
  raster::as.data.frame(xy = TRUE) %>% 
  as_tibble() %>%
  rename(lon = 1, lat = 2, chl = 3) %>%
  mutate(chl_log = log(chl + 1))  # Convert to log scale

## map the distribution of chl
ggplot() + 
  geom_raster(data = chl.a.tb, aes(x = lon, y = lat, fill = chl_log)) +
  geom_sf(data = world, fill = "grey80", col = "black") +
  coord_sf(xlim = c(-21, -13), c(24.2, 28.5)) +
  scale_fill_gradientn(colours = oce::oce.colorsJet(120), 
                       na.value = "white", 
                       breaks = seq(0, log(4.2 + 1), 0.5), 
                       labels = scales::number_format(scale = 1)) +
  guides(fill = guide_colorbar(title = expression(log(Chlorophyll_a + 1))),
         title.position = "right",
         title.theme = element_text(angle = 90), 
         barwidth = unit(.5, "cm"), barheight = unit(7.5, "cm"), title.hjust = .5) +
  theme_bw() 

THIS IS THE CURRENT OUTPUT:

THIS IS MY DESIRED OUTPUT: with higher resolution, log scale values and bar as legend.

Answer 1

here some hints using terra package (raster is now deprecated).

Border values (holes and boundaries) won't be interpolated because of NA's (in your desired plot it's assumed that they are zeroes!)

I did not put the rest of your code (for cropping and plot format) as it is basically the same.

• terra::disagg function is a simple resampling function for integer scale factor.
• tidyterra has convenient functions to help plot terra rasters (it does the same as as.data.frame(xy=TRUE))

library(terra)
library(tidyterra)
library(ggplot2)

# load original raster
test <- terra::rast("file_chl_satellite.nc")

# supersample by a factor of 4
test <- terra::disagg(test, 4, method="bilinear")

# clamp values
test <- terra::clamp(test, lower=0.01, upper=0.5, values=FALSE)

# rescale values
test <- log(test+1)

# rename layer (easier for setting fill aes)
names(test) <- "Chl"

# build plot
ggplot()+
  tidyterra::geom_spatraster(
    aes(fill=Chl), 
    data=test, 
    interpolate=TRUE # if you want further interpolation
  )+
  scale_fill_gradientn(
    colours = oce::oce.colorsJet(120),
    breaks = c(0.1, 0.2, 0.3, 0.4, 0.5),
    limits = c(0.1, 0.5)
  ) +
  theme(legend.position="bottom")+
  guides(
    fill = guide_colorbar(
      title = expression(log(Chlorophyll_a + 1)), 
      title.position = "bottom", 
      barwidth = unit(10, "cm"), 
      barheight = unit(.5, "cm"),
      title.hjust = .5,
      frame.colour="black",
      ticks.colour="black"
    )
  )