knitr::opts_chunk$set(echo = TRUE)
library(unittest) # Redirect ok() output to stderr ok <- function(...) capture.output(unittest::ok(...), file = stderr())
library(mfdb)
# Don't show chatty MFDB messages logging::setLevel('WARN')
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
# NB: Calculated with: # raw <- hydroVar[hydroVar$year == 2018 & hydroVar$rect %in% c('44G7','44G8','44G9','45G7','45G8','45G9') & hydroVar$month %in% c(1,2,3),] # sum(raw$sob * raw$areaSea_km2) / sum(raw$areaSea_km2) ok(ut_cmp_equal(dat, read.table(blank.lines.skip = TRUE, header = TRUE, stringsAsFactors = FALSE, text = ' year step area mean 2018 1 a1 5.131353 2018 1 a2 5.286023 2018 2 a1 5.772875 2018 2 a2 5.787779 2018 3 a1 7.259383 2018 3 a2 6.609905 2018 4 a1 6.144196 2018 4 a2 6.020029 ', colClasses = c(NA, 'character', NA, NA)), tolerance = 1e-5), "extract mean bottom temperature by quarter and area")
mfdb_disconnect(mdb)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.