Dsl_depth

Calc Z_pds**l from Proud

eqn s6, proud

Z_{pdsl*</sub> = 483.8 + (1272 * *tao*) − (143 * *PP)}

where:

Z_pds**l: depth of the principle DSL

P**P: primary production was converted from carbon to wet-mass biomass using the conversion factor 1:9 (g C m-2 day-1)

tao: wind stress (N m-2)

Get tao

#------------
## WIND -> TAO
# get wind data and calculate tao for input to Proud function for DSL

devtools::load_all('~/analyzePSAT')
library(curl)

limits <- list(xmin=-93.5, xmax=-4, ymin=-10, ymax=55)
tVec <- seq.Date(as.Date('2004-01-15'), Sys.Date(), by='month')
#time <- c(as.Date('2004-01-01'), Sys.Date())
#filename
#save.dir

for (i in 70:length(tVec)){
  filename <- paste('SeaWinds')
  save.dir <- '~/work/EnvData/wind/'
  time <- tVec[i]
  repeat{
    #get.hycom(spatLim, time, filename = paste(format(time, '%Y%m%d'), '_', filename, '.nc', sep = ''),
    #          download.file = TRUE, dir = save.dir, depLevels=depLevels, ...) 
    get.wind(limits=limits, 
             time=time, dir=save.dir,
             filename=paste(format(time, '%Y%m%d'), '_', filename, '.nc', sep = ''))
    tryCatch({
      err <- try(RNetCDF::open.nc(paste(save.dir,'/', format(time, '%Y%m%d'), '_', filename, '.nc', sep = '')), silent = T)
    }, error=function(e){print(paste('ERROR: Download of data at ', time, ' failed. Trying call to server again.', sep = ''))})
    if(class(err) != 'try-error') break
  }
}

Get VGPM

#------------
## VGPM PP
# see http://lukemiller.org/index.php/2011/12/loading-osus-vgpm-ocean-productivity-data-in-r/


setwd('~/work/EnvData/pp_vgpm/')
tVec <- seq.Date(as.Date('2004-01-01'), Sys.Date(), by='month')
for (i in 1:length(tVec)){
  doy <- lubridate::yday(tVec[i])
  if (nchar(doy) == 1) doy <- paste('00', doy, sep='')
  if (nchar(doy) == 2) doy <- paste('0', doy, sep='')
  yrday <- paste(lubridate::year(tVec[i]), doy, sep='')
  url <- paste('http://orca.science.oregonstate.edu/data/1x2/monthly/vgpm.r2018.m.chl.m.sst/xyz/vgpm.', yrday, '.all.xyz.gz', sep='')
  filename <- paste('vgpm.', yrday, '.all.xyz.gz', sep='')
  #repeat{
  curl::curl_download(url, filename, quiet=FALSE)
  
  #  tryCatch({
  #    err <- try(RNetCDF::open.nc(filename), silent = T)
  #  }, error=function(e){print(paste('ERROR: Download of data at ', yrday, ' failed. Trying call to server again.', sep = ''))})
  #  if(class(err) != 'try-error') break
  #}
  
}

Calculate Z_pds**l

limits <- list(xmin=-93.5, xmax=-4, ymin=-10, ymax=55)
tVec.wind <- seq.Date(as.Date('2004-01-15'), Sys.Date(), by='month')
tVec.pp <- seq.Date(as.Date('2004-01-01'), Sys.Date(), by='month')

for (i in 20:length(tVec.pp)){
  
  # get proper date formating for vgpm
  doy <- lubridate::yday(tVec.pp[i])
  if (nchar(doy) == 1) doy <- paste('00', doy, sep='')
  if (nchar(doy) == 2) doy <- paste('0', doy, sep='')
  yrday <- paste(lubridate::year(tVec.pp[i]), doy, sep='')
  
  # load vgpm and rasterize
  t1 <- Sys.time()
  vgpm <- vgpm.load(file=paste('~/work/EnvData/pp_vgpm/vgpm.', yrday, '.all.xyz', sep=''), w.lon = -180, e.lon = 180,
                    n.lat = 90, s.lat = -90, raster = TRUE)
  vgpm <- vgpm / 1000 # convert from mg to g C/ m^2 / d^1
  t2 <- Sys.time()
  print(paste(t2 - t1)) # print how long it took
  
  # load wind raster. units are 1 Pa = 1 N / m2
  taux <- raster(paste('~/work/EnvData/wind/', format(tVec.wind[i], '%Y%m%d'), '_SeaWinds.nc', sep = ''), varname='taux')
  tauy <- raster(paste('~/work/EnvData/wind/', format(tVec.wind[i], '%Y%m%d'), '_SeaWinds.nc', sep = ''), varname='tauy')
  tau <- sqrt(taux^2 + tauy^2)
  rm(taux); rm(tauy)
  
  # calculate Z_pdsl
  vgpm <- resample(vgpm, tau)
  vgpm <- vgpm / cellStats(vgpm, 'max')
  zpdsl <- 483.8 + 1272 * tau - 143 * vgpm
  zpdsl[zpdsl < 0] <- 0
  writeRaster(zpdsl, paste('~/work/EnvData/zpdsl/zpdsl_', tVec.pp[i], '.grd', sep=''))
  
  #rm(tau); rm(vgpm); gc()
}

setwd('~/work/RCode/deepdiving/dsl_depth_notebook/')
fList <- list.files()
fList <- fList[grep('zpdsl', fList)]
fList <- fList[grep('grd', fList)]
#setwd('./dsl_depth_notebook/')
s <- stack(fList)
names(s) <- c('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec')

plot(s, zlim=c(400,900), main='Monthly DSL depth in 2004 following Proud')

Calc Z_dvm from Bianchi

bianchi supp eqn

Z_DVM = 398 − 0.56 * deltaO* − 115 * *log*10(*chl*) + 0.36 * *mld* − 2.4 * *deltaT

where:

delta**O: oxygen diff between 0-25m and 150-500m (mmol m-3)

delta**T: temp diff between 0-25m and 150-500m (deg C)

chl: log transformed surface CHL (mg m-3)

mld: mixed layer depth (m). deltaT=0.2degC from 10m (not surface)

Monthly CHL

#------------
## CHL

setwd('~/work/EnvData/chl/')
for (bb in 2004:2018){
  url <- paste('ftp://my.cmems-du.eu/Core/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/dataset-oc-glo-chl-multi-l4-gsm_25km_monthly-rep-v02/', bb, '/', sep='')
  fList <- unlist(strsplit(RCurl::getURL(url, verbose=TRUE, ftp.use.epsv=TRUE, dirlistonly = TRUE, userpwd = 'cbraun:CamrinCMEMS2017'), '\n'))
  for (i in 1:length(fList)){
    url <- paste('ftp://cbraun:CamrinCMEMS2017@my.cmems-du.eu/Core/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/dataset-oc-glo-chl-multi-l4-gsm_25km_monthly-rep-v02/', bb, '/', sep='')
    getFile <- paste(url, fList[i], sep='')
    download.file(getFile, fList[i])
  }
  
}

Monthly oxygen at depth

#------------
## OXYGEN
# comes from the motu python client but cant run using system() for some reason. download
# extent as .nc using python then extract from raster stack here
for (i in 1:length(allList)){
  if(i==1){
    all <- allList[[i]]$track
  } else{
    all <- rbind(all, allList[[i]]$track)
  }
}
all$xdist[which(all$xdist >=2)] <- 2
all$ydist[which(all$ydist >=2)] <- 2
all$xdist[which(is.na(all$xdist))] <- 1
all$ydist[which(is.na(all$ydist))] <- 1

## consider extracting oxygen by month
## this commented one works for sure
#python ~/Documents/motu-client-python/motu-client.py --user cbraun --pwd CamrinCMEMS2017 --motu http://my.cmems-du.eu/motu-web/Motu --service-id GLOBAL_REANALYSIS_BIO_001_018-TDS --product-id dataset-global-nahindcast-bio-001-018-V5-o2 --longitude-min -80 --longitude-max -75 --latitude-min 30 --latitude-max 35 --date-min "2016-12-16 12:00:00" --date-max "2016-12-16 12:00:00" --depth-min 0.5056 --depth-max 0.5059 --variable O2 --variable nav_lon --variable nav_lat --out-dir ~/Documents/ --out-name try-motu.nc
system('python ~/Documents/motu-client-python/motu-client.py --user cbraun --pwd CamrinCMEMS2017 --motu http://my.cmems-du.eu/motu-web/Motu --service-id GLOBAL_REANALYSIS_BIO_001_018-TDS --product-id dataset-global-nahindcast-bio-001-018-V5-o2 --longitude-min -94 --longitude-max -4 --latitude-min -10 --latitude-max 49 --date-min 2004-09-24 00:00:00 --date-max 2009-09-31 00:00:00 --depth-min 0.5056 --depth-max 508.64 --variable O2 --variable nav_lon --variable nav_lat --out-dir ~/Documents/ --out-name motu-o2-092004.nc')
#pycall <- paste('python ~/Documents/motu-client-python/motu-client.py --user cbraun --pwd CamrinCMEMS2017 --motu http://my.cmems-du.eu/motu-web/Motu --service-id GLOBAL_REANALYSIS_BIO_001_018-TDS --product-id dataset-global-nahindcast-bio-001-018-V5-o2 --longitude-min ', xlims[1],
#' --longitude-max ', xlims[2], ' --latitude-min 30 --latitude-max 35 --date-min "2016-12-16 12:00:00" --date-max "2016-12-16 12:00:00" --depth-min 0.5056 --depth-max 0.5059 --variable O2 --variable nav_lon --variable nav_lat --out-dir ~/Documents/ --out-name try-motu.nc)

# eventually need deltaO
limits <- list(xmin=-93.5, xmax=-4, ymin=-10, ymax=55)
tVec.o2 <- seq.POSIXt(as.POSIXct('2004-01-16 12:00:00'), as.POSIXct('2018-01-16 12:00:00'), by='month')

for (i in 98:length(tVec.o2)){
  system(paste('python ~/Documents/motu-client-python/motu-client.py --user cbraun --pwd CamrinCMEMS2017 --motu http://my.cmems-du.eu/motu-web/Motu --service-id GLOBAL_REANALYSIS_BIO_001_029-TDS --product-id global-reanalysis-bio-001-029-monthly --longitude-min ', limits$xmin,' --longitude-max ', limits$xmax, ' --latitude-min ', limits$ymin, ' --latitude-max ', limits$ymax, ' --date-min "', format(tVec.o2[i]-36*60*60-1, '%Y-%m-%d %H:%M:%S'), '" --date-max "', format(tVec.o2[i]+36*60*60+1, '%Y-%m-%d %H:%M:%S'), '" --depth-min 0.5056 --depth-max 628.0269999999999 --variable o2 --out-dir ~/work/EnvData/oxygen --out-name ', format(tVec.o2[i], '%Y%m%d'),'_o2.nc', sep=''))
}

fList <- list.files('~/work/EnvData/oxygen/')
fList <- substr(fList, 1, 6)
tVec.o2.mod <- tVec.o2[-which(format(tVec.o2, '%Y%m') %in% fList)]

for (i in 1:length(tVec.o2.mod)){
  system(paste('python ~/Documents/motu-client-python/motu-client.py --user cbraun --pwd CamrinCMEMS2017 --motu http://my.cmems-du.eu/motu-web/Motu --service-id GLOBAL_REANALYSIS_BIO_001_029-TDS --product-id global-reanalysis-bio-001-029-monthly --longitude-min ', limits$xmin,' --longitude-max ', limits$xmax, ' --latitude-min ', limits$ymin, ' --latitude-max ', limits$ymax, ' --date-min "', format(tVec.o2.mod[i]-120*60*60-1, '%Y-%m-%d %H:%M:%S'), '" --date-max "', format(tVec.o2.mod[i]+120*60*60+1, '%Y-%m-%d %H:%M:%S'), '" --depth-min 0.5056 --depth-max 628.0269999999999 --variable o2 --out-dir ~/work/EnvData/oxygen --out-name ', format(tVec.o2.mod[i], '%Y%m%d'),'_o2.nc', sep=''))
}

fList <- list.files('~/work/EnvData/oxygen/')
fList <- substr(fList, 1, 6)
tVec.o2[-which(format(tVec.o2, '%Y%m') %in% fList)]
#tVec.o2.mod <- 

url <- 'ftp://my.cmems-du.eu/Core/OCEANCOLOUR_GLO_CHL_L4_REP_OBSERVATIONS_009_082/dataset-oc-glo-chl-multi-l4-gsm_25km_monthly-rep-v02/2015/'#20150101_m_20150131-ACRI-L4-CHL-MULTI_25KM-GLO-REP-v02.nc
fList <- unlist(strsplit(RCurl::getURL(url, verbose=TRUE, ftp.use.epsv=TRUE, dirlistonly = TRUE, userpwd = 'cbraun:CamrinCMEMS2017'), '\n'))
for (i in 1:length(fList)){
  download.file()
}

Monthly temp at depth

#------------
## MONTHLY TEMP AT DEPTH FOR BIANCHI

library(RCMEMS)

setwd('~/work/EnvData/tz_proud')

limits <- list(xmin=-93.5, xmax=-4, ymin=-10, ymax=55)
tVec.monthTZ <- seq.POSIXt(as.POSIXct('2004-01-16 12:00:00'), as.POSIXct('2018-08-16 12:00:00'), by='month')
for (i in 2:length(tVec.monthTZ)){
  #system(paste('python ~/Documents/motu-client-python/motu-client.py --user cbraun --pwd CamrinCMEMS2017 --motu http://my.cmems-du.eu/motu-web/Motu --service-id GLOBAL_REANALYSIS_BIO_001_029-TDS --product-id global-reanalysis-bio-001-029-monthly --longitude-min ', limits$xmin,' --longitude-max ', limits$xmax, ' --latitude-min ', limits$ymin, ' --latitude-max ', limits$ymax, ' --date-min "', format(tVec.monthTZ[i]-120*60*60-1, '%Y-%m-%d %H:%M:%S'), '" --date-max "', format(tVec.monthTZ[i]+120*60*60+1, '%Y-%m-%d %H:%M:%S'), '" --depth-min 0.5056 --depth-max 628.0269999999999 --variable o2 --out-dir ~/work/EnvData/oxygen --out-name ', format(tVec.monthTZ[i], '%Y%m%d'),'_o2.nc', sep=''))
  
  cfg <- CMEMS.config(python='/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python',
                      script='~/Documents/motuclient-python/motuclient.py',
                      user='cbraun',
                      pwd='CamrinCMEMS2017',
                      motu='http://my.cmems-du.eu/motu-web/Motu',
                      service.id = 'GLOBAL_REANALYSIS_PHY_001_030-TDS',
                      product.id = 'global-reanalysis-phy-001-030-monthly',
                      variable = 'thetao',
                      out.dir = getwd(),
                      out.name = paste(format(tVec.monthTZ[i], '%Y%m%d'),'_mTZ.nc', sep=''))
  
  CMEMS.download(cfg,
                 date.min = format(tVec.monthTZ[i]-120*60*60+1, '%Y-%m-%d %H:%M:%S'),
                 date.max = format(tVec.monthTZ[i]+120*60*60+1, '%Y-%m-%d %H:%M:%S'),
                 latitude.min = limits$ymin,
                 latitude.max = limits$ymax,
                 longitude.min = limits$xmin,
                 longitude.max = limits$xmax,
                 depth.min = '0',
                 depth.max = '550',
                 debug=FALSE)
  
  cfg <- CMEMS.config(python='/Library/Frameworks/Python.framework/Versions/2.7/Resources/Python.app/Contents/MacOS/Python',
                      script='~/Documents/motuclient-python/motuclient.py',
                      user='cbraun',
                      pwd='CamrinCMEMS2017',
                      motu='http://my.cmems-du.eu/motu-web/Motu',
                      service.id = 'GLOBAL_REANALYSIS_PHY_001_030-TDS',
                      product.id = 'global-reanalysis-phy-001-030-monthly',
                      variable = 'mlotst',
                      out.dir = getwd(),
                      out.name = paste(format(tVec.monthTZ[i], '%Y%m%d'),'_mld.nc', sep=''))
  
  CMEMS.download(cfg,
                 date.min = format(tVec.monthTZ[i]-120*60*60+1, '%Y-%m-%d %H:%M:%S'),
                 date.max = format(tVec.monthTZ[i]+120*60*60+1, '%Y-%m-%d %H:%M:%S'),
                 latitude.min = limits$ymin,
                 latitude.max = limits$ymax,
                 longitude.min = limits$xmin,
                 longitude.max = limits$xmax,
                 depth.min = '0',
                 depth.max = '550',
                 debug=FALSE)
  
  
  
}

Calc Z_dvm from Bianchi

#------------
## CALC DVM FROM BIANCHI

tVec <- seq.POSIXt(as.POSIXct('2004-01-01'), as.POSIXct('2016-12-01'), by='month')

for (i in 1:length(tVec)){
  #for (i in 1:3){
  # get chl
  fList <- list.files('~/work/EnvData/chl/dvm_calc/', full.names = T)
  fList <- fList[grep(format(tVec[i], '%Y%m'), fList)]
  chl <- raster(fList)
  
  # get o2
  fList <- list.files('~/work/EnvData/oxygen/', full.names = T)
  fList <- fList[grep(format(tVec[i], '%Y%m'), fList)]
  o2 <- brick(fList, lvar=4)
  if (i == 1){
    nc <- RNetCDF::open.nc(fList[1])
    depth_o2 <- RNetCDF::var.get.nc(nc, 'depth')
  }
  #upper_o2 <- mean(o2[[1:14]], na.rm=T) # 0-25 m
  #lower_o2 <- mean(o2[[28:40]], na.rm=T) # 150-500 m
  delta_o2 <- mean(o2[[1:14]], na.rm=T) - mean(o2[[28:40]], na.rm=T)
  
  # get T
  fList <- list.files('~/work/EnvData/monthlyTZ/', full.names = T)
  fList <- fList[grep('mTZ', fList)]
  fList <- fList[grep(format(tVec[i], '%Y%m'), fList)]
  temps <- brick(fList, lvar=4)
  if (i == 1){
    nc <- RNetCDF::open.nc(fList)
    depth_t <- RNetCDF::var.get.nc(nc, 'depth')
    
  }
  #upper_t <- mean(temps[[1:14]], na.rm=T) # 0-25 m
  #lower_t <- mean(temps[[25:32]], na.rm=T) # 150-500 m
  delta_t <- mean(temps[[1:14]], na.rm=T) - mean(temps[[25:32]], na.rm=T)
  
  # get MLD
  fList <- list.files('~/work/EnvData/monthlyTZ/', full.names = T)
  fList <- fList[grep('mld', fList)]
  fList <- fList[grep(format(tVec[i], '%Y%m'), fList)]
  mld <- raster(fList)
  
  # resample everything to most coarse (0.25) and cut to study area (n atl)
  mld <- raster::resample(mld, delta_o2)
  chl <- raster::resample(chl, delta_o2)
  delta_t <- raster::resample(delta_t, delta_o2)
  
  z_dvm <- 398 - 0.56 * delta_o2 - 115 * log10(chl) + 0.36 * mld - 2.4 * delta_t
  
  #Z_DVM = 398 - 0.56 * deltaO2 - 115* Log10(Chl) + 0.36 * mld - 2.4 * deltaT
  #where deltaO2 = data.o2grad:
  #  data.o2grad = double(squeeze(nanmean(data.o2(1:4,:,:),1) - nanmean(data.o2(9:14,:,:),1)));
  # deltao2 <- drop(mean(data.o2[1:4,,], na.rm=T) - mean(data.o2[9:14,,], na.rm=T))
  
  setwd('~/work/EnvData/zdvm/')
  writeRaster(z_dvm, filename=paste(format(tVec[i], '%Y%m'), '_zdvm.grd',sep=''), overwrite=T)
  
}

Perform just an example calc

## get monthly data for 2004 as an example
#system('aws s3 cp s3://braun-data/EnvData/chl/dvm_calc/ . --exclude "*" --include "*2004*" --recursive')
#system('aws s3 cp s3://braun-data/EnvData/oxygen/ . --exclude "*" --include "*2004*" --recursive')
#system('aws s3 cp s3://braun-data/EnvData/monthlyTZ/ . --exclude "*" --include "*2004*" --recursive')

fList <- list.files()

fList_chl <- fList[grep('CHL', fList)]
fList_o2 <- fList[grep('o2', fList)]
fList_temp <- fList[grep('mTZ', fList)]
fList_mld <- fList[grep('mld', fList)]

for (i in 1:7){
  chl <- raster(fList_chl[i]) #chl/dvm_calc
  
  o2 <- brick(fList_o2[i], lvar=4)
  delta_o2 <- mean(o2[[1:14]], na.rm=T) - mean(o2[[28:40]], na.rm=T)
  
  temps <- brick(fList_temp[i], lvar=4)
  delta_t <- mean(temps[[1:14]], na.rm=T) - mean(temps[[25:32]], na.rm=T)
  
  mld <- raster(fList_mld[i])
  
  # resample everything to most coarse (0.25) and cut to study area (n atl)
  mld <- raster::resample(mld, delta_o2)
  chl <- raster::resample(chl, delta_o2)
  delta_t <- raster::resample(delta_t, delta_o2)

  if (i == 1){
      z_dvm <- 398 - 0.56 * delta_o2 - 115 * log10(chl) + 0.36 * mld - 2.4 * delta_t
      names(z_dvm) <- 'jan'
  } else{
      z_dvm <- stack(z_dvm, 398 - 0.56 * delta_o2 - 115 * log10(chl) + 0.36 * mld - 2.4 * delta_t)
      mnames <- c('jan','feb','mar','apr','may','june','july','aug')
      names(z_dvm)[i] <- mnames[i]
  }
}

plot(z_dvm, zlim=c(200,600), main='Monthly DSL depth in 2004 following Bianchi')

Look at the difference for select months

proud <- raster::resample(s[[1:7]], z_dvm)
diffs <- proud - z_dvm
names(diffs) <- c('jan','feb','mar','apr','may','jun','jul')
plot(diffs, zlim=c(0, 600), main='Proud minus Bianchi for Jan thru July 2004')

Calc Zpds**l from Proud

Get tao

Get VGPM

Calculate Zpds**l

Calc Zdvm from Bianchi

Monthly CHL

Monthly oxygen at depth

Monthly temp at depth

Calc Zdvm from Bianchi

Perform just an example calc

Look at the difference for select months

Calc Z_pds**l from Proud

Calculate Z_pds**l

Calc Z_dvm from Bianchi

Calc Z_dvm from Bianchi