library(mfdb)Hydrographic products are available from numerous sources. Here, we provide a working example where outputs from a hydrographic model are:
CMEMS Baltic Sea Physical Reanalysis product provides a physical reanalysis for the whole Baltic Sea area from January 1993 and up to minus 1-1.5 year compared to real time. The product is produced by using the ice-ocean model NEMO-Nordic (based on NEMO-3.6, Nucleus for European Modelling of the Ocean) together with a LSEIK data assimilation scheme, https://resources.marine.copernicus.eu/?option=com_csw&view=details&product_id=BALTICSEA_REANALYSIS_PHY_003_011. As a demonstration we retrieve one year of already processed monthly means of bottom temperature (bottomT), but hydrographic products are also available for separate depth layers and at higher spatiotemporal resolution, and can be aggregated before import according to the needs. Note that an account has to be created to download data from Copernicus Marine Service using the script below.
library(RCMEMS)
# point to the data product
cfg <- CMEMS.config(motu="http://my.cmems-du.eu/motu-web/Motu",
service.id = "BALTICSEA_REANALYSIS_PHY_003_011-TDS",
product.id = "dataset-reanalysis-nemo-monthlymeans",
variable = c("bottomT"))
# add user and psw
CMEMS.config.usr <- function(x){
print("username")
scan("", what="character", nmax=1, quiet=T)
}
CMEMS.config.pwd <- function(x){
print("password")
scan("", what="character", nmax=1, quiet=T)
}
cfg@user <- CMEMS.config.usr()
cfg@pwd <- CMEMS.config.pwd()
# select year and download
y <- 2018 # year of interest
CMEMS.download(cfg,
ROI = c(17,20,56,58.5),
date.range = c(ISOdate(y,01,01), ISOdate(y,12,31)),
depth.range= c(10,500), # max depth Baltic 459 m
out.path=paste("data_provided/CMEMS_BAL_PHY_reanalysis_monthlymeans_",y,".nc",sep=""),
debug=FALSE)temp <- tempfile()
download.file(url="https://gis.ices.dk/shapefiles/ICES_StatRec_mapto_ICES_Areas.zip", temp)
unzip(temp, exdir="data_provided/ices_rect")
unlink(temp)library(sf)
library(raster)
library(tidyverse)
library(rnaturalearth)
library(lwgeom)
library(ggplot2)
y <- 2018 # year of interest
# load ICES rect
ices.rect <- st_read("data_provided/ices_rect/StatRec_map_Areas_Full_20170124.dbf", quiet = T)
# load bottom temperature (bottomT)
rst <- brick(paste("data_provided/CMEMS_BAL_PHY_reanalysis_monthlymeans_",y,".nc",sep=""), varname="bottomT", lvar=4)
# calculate mean sob (raster) by ices rect (polygons)
ov <- raster::extract(rst,ices.rect, fun=mean, na.rm=T, df=T)
hydroVar <- data.frame(ID=1:length(ices.rect$ICESNAME),
rect=ices.rect$ICESNAME,
areaFull_km2=ices.rect$AREA_KM2,
year=y) %>%
right_join(ov) %>%
select(-ID) %>%
gather("month","sob",4:15) %>%
mutate(month = as.numeric(substring(month,7,8))) %>%
subset(!is.na(sob)) # NA are on land or outside the raster area
land <- rnaturalearth::ne_countries(returnclass = "sf", continent="Europe", scale="large") %>%
st_union()
ices.rect.sea <- ices.rect %>% filter(ICESNAME %in% hydroVar$rect) %>%
st_difference(land)
ices.rect.sea$areaSea_km2 <- st_area(ices.rect.sea) %>% units::set_units(km^2)
hydroVar <- ices.rect.sea %>%
rename(rect=ICESNAME) %>%
select(rect, areaSea_km2) %>%
right_join(hydroVar)
ggplot(hydroVar) +
geom_sf() +
geom_sf(data=st_geometry(ices.rect.sea), fill="lightblue") +
xlim(xmin(extent(hydroVar)), xmax(extent(hydroVar))) +
ylim(ymin(extent(hydroVar)), ymax(extent(hydroVar))) +
theme_bw()
# Write the results to temporary files so we can read them back in the next step
write.table(sf::st_drop_geometry(ices.rect.sea), file = 'ices.rect.sea.txt')
write.table(sf::st_drop_geometry(hydroVar), file = 'hydroVar.txt')Import ICES rectangles and import the hydrographic data as a survey index. Then we can retrieve the hydrographic data as a weighted mean using ICES rectangle as a weighting variable:
mdb <- mfdb(tempfile(fileext = '.duckdb'))
ices.rect.sea <- read.table('ices.rect.sea.txt')
hydroVar <- read.table('hydroVar.txt')
# import ICES rectangles
mfdb_import_area(mdb, data.frame(
name = as.character(ices.rect.sea$ICESNAME),
size = ices.rect.sea$areaSea_km2))
# import area size as index for later use
mfdb_import_cs_taxonomy(mdb, "index_type", data.frame(name=c("ices_rect", "bottom_temp")))
mfdb_import_survey_index(mdb,
data_source="area_ices_rect",
data.frame(index_type = "ices_rect",
year = hydroVar$year,
month = hydroVar$month,
areacell = as.character(hydroVar$rect),
value = as.numeric(hydroVar$areaSea_km2)))
# import avg bottom temperature by month and ICES rect (under the 'length' field)
mfdb_import_survey_index(mdb,
data_source="bottom_temp_CMEMS",
data.frame(year = hydroVar$year,
month = hydroVar$month,
areacell = as.character(hydroVar$rect),
value = hydroVar$sob,
index_type = "bottom_temp",
stringsAsFactors = TRUE))
# extract mean bottom temperature by quarter and area (weighted mean using rectangle area)
dat <- mfdb_survey_index_mean(mdb, c(), list(
## area = mfdb_unaggregated(),
area = mfdb_group('a1'=c('42G7','42G8','42G9','43G7','43G8','43G9'),
'a2'=c('44G7','44G8','44G9','45G7','45G8','45G9')),
timestep = mfdb_timestep_quarterly,
year = 2018,
index_type = "bottom_temp",
data_source = 'bottom_temp_CMEMS'),
scale_index = 'area_size')[[1]]
dat <- dat[,c('year', 'step', 'area', 'mean')]
dat## year step area mean
## 1 2018 1 a1 5.131353
## 2 2018 1 a2 5.286023
## 3 2018 2 a1 5.772875
## 4 2018 2 a2 5.787779
## 5 2018 3 a1 7.259383
## 6 2018 3 a2 6.609905
## 7 2018 4 a1 6.144196
## 8 2018 4 a2 6.020029mfdb_disconnect(mdb)