rggplot2rastersatellitencdf4

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


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:

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: Plot with too many pixels

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


Solution

  • 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.

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