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 MY DESIRED OUTPUT: with higher resolution, log scale values and bar as legend.
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"
)
)