rmd/reprojectAMM7.R

In NERC-CEH/ukem: Package Title (one line, title case)

# This script reads the ocean model AMM7 output for 
# CO2 and N2O fluxes (netCDF files).
# This is cropped to the NAME grid, daily output is 
# averaged to annual and/or monthly means, and 
# resampled on to the NAME grid.

# cd /gws/nopw/j04/dare_uk/plevy/AMM7/
# module load jaspy/3.7/r20200606
rm(list=ls(all=TRUE))
library(raster)
library(spCEH)
library(ggplot2)
library(lubridate)
library(ncdf4)
library(terra)

"%!in%" <- Negate("%in%")
Sys.setenv(TZ = "UTC")
rasterOptions(maxmemory=1e8, datatype = "FLT4S")
rasterOptions(tmpdir = "/work/scratch-nopw/bkruijt/R")
r_name <- getRasterTemplate("UK_NAME", proj = projlonlat)

beginCluster(n = 12)
# setwd("E:/plevy/AMM7")
# setwd("C:/Users/plevy/OneDrive - NERC/0Peter/curr/DARE/AMM7")

cropAvgResample <- function(fname, varname, r, averaging = "annual", resampling = TRUE){
  b <- brick(fname, varname = varname)
  # get dimensions
  ncfile <- nc_open(fname)
  time <- ncvar_get(ncfile, "time_counter")
  lon <- ncvar_get(ncfile, "lon")
  lat <- ncvar_get(ncfile, "lat")
  nc_close(ncfile)

  # get datetime vector
  datect <- as.POSIXct(time, origin = "1950-01-01")
  datelt <- as.POSIXlt(datect)

  # get spatial grid dimesnsions: resolution and extent
  res = 0.11111
  ext <- extent(min(lon) - res/2, max(lon) + res/2,
                min(lat) - res/2, max(lat) + res/2)
  extent(b) <- ext
  crs(b) <- projlonlat
  
  b <- crop(b, r)
  
  # calc annual or monthly averages
  if (averaging == "annual"){
    b <- stackApply(b, indices=datelt$year-104, fun=mean)
    names(b) <- 2005:2015
  }
  if (averaging == "monthly"){
    b <- stackApply(b, indices=datelt$mon+1, fun=mean)
    names(b) <- 1:12
  }
  
  # resample to new raster grid
  if (resampling) b <- terra::resample(b, r)
  return(b)
}

### co2 ###
# some confusion over units: either mmol / m2 / d or mol / m2 / s
#  mult <- 1/1000/86400 # but has this been applied already or not?
b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-CO2flux_2005-2015.nc",
  varname = "CO2_sea_air_flux", r = r_name, averaging = "annual")
# write to tif file
rf <-  writeRaster(b, filename = "b_Fco2_annual.tif", overwrite=TRUE) 
# get overall mean
b <- mean(b)
rf <-  writeRaster(b, filename = "b_Fco2_mean.tif", overwrite=TRUE) 

b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-CO2flux_2005-2015.nc",
    varname = "CO2_sea_air_flux", r = r_name, averaging = "monthly")
rf <-  writeRaster(b, filename = "b_Fco2_monthly.tif", overwrite=TRUE) 

### n2o ###
b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-N2Oflux_2005-2015.nc",
  varname = "N2O_sea_air_flux", r = r_name, averaging = "annual")
# write to tif file
rf <-  writeRaster(b, filename = "b_Fn2o_annual.tif", overwrite=TRUE) 
# get overall mean
b <- mean(b)
rf <-  writeRaster(b, filename = "b_Fn2o_mean.tif", overwrite=TRUE) 

b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-N2Oflux_2005-2015.nc",
    varname = "N2O_sea_air_flux", r = r_name, averaging = "monthly")
rf <-  writeRaster(b, filename = "b_Fn2o_monthly.tif", overwrite=TRUE) 


###  repeat but without resampling ###
### co2 ###
b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-CO2flux_2005-2015.nc",
  varname = "CO2_sea_air_flux", r = r_name, averaging = "annual", resampling = FALSE)
rf <-  writeRaster(b, filename = "b_Fco2_annual_amm7.tif", overwrite=TRUE) 

b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-CO2flux_2005-2015.nc",
    varname = "CO2_sea_air_flux", r = r_name, averaging = "monthly", resampling = FALSE)
rf <-  writeRaster(b, filename = "b_Fco2_monthly_amm7.tif", overwrite=TRUE) 

### n2o ###
b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-N2Oflux_2005-2015.nc",
  varname = "N2O_sea_air_flux", r = r_name, averaging = "annual", resampling = FALSE)
rf <-  writeRaster(b, filename = "b_Fn2o_annual_amm7.tif", overwrite=TRUE) 

b <- cropAvgResample(fname <- "../../AMM7/DAREUK_AMM7-N2Oflux_2005-2015.nc",
    varname = "N2O_sea_air_flux", r = r_name, averaging = "monthly", resampling = FALSE)
rf <-  writeRaster(b, filename = "b_Fn2o_monthly_amm7.tif", overwrite=TRUE) 


endCluster()