R/temperature.db.r

temperature.db = function ( ip=NULL, p, DS, varnames=NULL, yr=NULL, ret="NULL", voi="t" ) {

  # over-ride default dependent variable name if it exists
  if (exists("variables",p)) if(exists("Y", p$variables)) voi=p$variables$Y
  
  if (DS=="lbm.inputs") {

    B = hydro.db( p=p, DS="bottom.all"  ) 
    B = B[ which(B$yr %in% p$yrs), ]
    B$tiyr = lubridate::decimal_date ( B$date )

    # globally remove all unrealistic data
    keep = which( B$t >= -3 & B$t <= 25 ) # hard limits
    if (length(keep) > 0 ) B = B[ keep, ]
    TR = quantile(B$t, probs=c(0.0005, 0.9995), na.rm=TRUE ) # this was -1.7, 21.8 in 2015
    keep = which( B$t >=  TR[1] & B$t <=  TR[2] )
    if (length(keep) > 0 ) B = B[ keep, ]
    keep = which( B$z >=  1 ) # ignore very shallow areas ..
    if (length(keep) > 0 ) B = B[ keep, ]
  
    varstokeep = unique( c( p$variables$Y, p$variables$LOCS, p$variables$TIME, p$variables$COV ) )
    B = B[,varstokeep]
    
    # default output grid
    Bout = bathymetry.db( p, DS="baseline", varnames=p$varnames )
    coords = p$variables$LOCS
    covars = p$variables$COV 
    if (length(covars)==1) {
      covs = list( Bout[,covars] )
      names(covs) = covars
      OUT  = list( LOCS = Bout[,coords], COV=covs ) 
    } else {
      OUT  = list( LOCS = Bout[,coords], COV=as.list( Bout[,covars] ) ) 
    }

    return (list(input=B, output=OUT))
  }

  # -----------------


  if ( DS %in% c("predictions", "predictions.redo" ) ) {
    # NOTE: the primary interpolated data were already created by lbm. 
    # This routine points to this data and also creates 
    # subsets of the data where required, determined by "spatial.domain.subareas" 
 
    outdir = file.path(project.datadirectory("bio.temperature"), "modelled", voi, p$spatial.domain )
    
    if (DS %in% c("predictions")) {
      P = V = NULL
      fn = file.path( outdir, paste("lbm.prediction", ret,  yr, "rdata", sep=".") )
      if (is.null(ret)) ret="mean"
      if (file.exists(fn) ) load(fn) 
      if (ret=="mean") return (P)
      if (ret=="sd") return( V)
    }

    if (exists( "libs", p)) RLibrary( p$libs )
    if ( is.null(ip) ) ip = 1:p$nruns

    # downscale and warp from p(0) -> p1

    for ( r in ip ) {
      # print (r)
      yr = p$runs[r, "yrs"]
      # default domain
      PP0 = lbm_db( p=p, DS="lbm.prediction", yr=yr, ret="mean")
      VV0 = lbm_db( p=p, DS="lbm.prediction", yr=yr, ret="sd")
      p0 = spatial_parameters( p=p, type=p$spatial.domain ) # from
      L0 = bathymetry.db( p=p0, DS="baseline" )
      L0i = lbm::array_map( "xy->2", L0[, c("plon", "plat")], gridparams=p0$gridparams )
      sreg = setdiff( p$spatial.domain.subareas, p$spatial.domain ) 

      for ( gr in sreg ) {
        p1 = spatial_parameters( p=p, type=gr ) # 'warping' from p -> p1
        L1 = bathymetry.db( p=p1, DS="baseline" )
        L1i = lbm::array_map( "xy->2", L1[, c("plon", "plat")], gridparams=p1$gridparams )
        L1 = planar2lonlat( L1, proj.type=p1$internal.crs )
        L1$plon_1 = L1$plon # store original coords
        L1$plat_1 = L1$plat
        L1 = lonlat2planar( L1, proj.type=p0$internal.crs )
        p1$wght = fields::setup.image.smooth( nrow=p1$nplons, ncol=p1$nplats, dx=p1$pres, dy=p1$pres, theta=p1$pres, xwidth=4*p1$pres, ywidth=4*p1$pres )

        P = matrix( NA, ncol=p$nw, nrow=nrow(L1) )
        V = matrix( NA, ncol=p$nw, nrow=nrow(L1) )
        for (iw in 1:p$nw) {
          P[,iw] = spatial_warp( PP0[,iw], L0, L1, p0, p1, "fast", L0i, L1i )
          V[,iw] = spatial_warp( VV0[,iw], L0, L1, p0, p1, "fast", L0i, L1i )
        }
        outdir_p1 = file.path(project.datadirectory("bio.temperature"), "modelled", voi, p1$spatial.domain)
        dir.create( outdir_p1, recursive=T, showWarnings=F )
        fn1_sg = file.path( outdir_p1, paste("lbm.prediction.mean",  yr, "rdata", sep=".") )
        fn2_sg = file.path( outdir_p1, paste("lbm.prediction.sd",  yr, "rdata", sep=".") )
        save( P, file=fn1_sg, compress=T )
        save( V, file=fn2_sg, compress=T )
        print (fn1_sg)

      }
    } 
    return ("Completed")

    if (0) {
      aoi = which( PS$z > 5 & PS$z < 3000 & PS$z.range < 500)
      levelplot( log(z) ~ plon + plat, PS[ aoi,], aspect="iso", labels=FALSE, pretty=TRUE, xlab=NULL,ylab=NULL,scales=list(draw=FALSE) )
      levelplot( log(t.ar_1) ~ plon + plat, PS[ aoi,], aspect="iso", labels=FALSE, pretty=TRUE, xlab=NULL,ylab=NULL,scales=list(draw=FALSE) )
      
      levelplot( log(t.range) ~ plon + plat, PS[ aoi,], aspect="iso", labels=FALSE, pretty=TRUE, xlab=NULL,ylab=NULL,scales=list(draw=FALSE) )
      levelplot( Z.rangeSD ~ plon + plat, PS[aoi,], aspect="iso", labels=FALSE, pretty=TRUE, xlab=NULL,ylab=NULL,scales=list(draw=FALSE) )
    }
  
  }


  # -----------------


  if (DS %in% c(  "bottom.statistics.annual", "bottom.statistics.annual.redo", "bottom.statistics.climatology" )){

		tstatdir = file.path( project.datadirectory("bio.temperature"), "modelled", voi, p$spatial.domain )
    dir.create( tstatdir, showWarnings=F, recursive = TRUE )

    if (DS == "bottom.statistics.climatology" ) {
      clim = NULL
      fn.climatology = file.path( tstatdir, paste("bottom.statistics.climatology", p$spatial.domain, "rdata", sep=".") )
      if (file.exists( fn.climatology ) ) load(fn.climatology)
      return ( clim )
    }

		if (DS %in% c("bottom.statistics.annual")) {
      BS = NULL
      fn = file.path( tstatdir, paste("bottom.statistics.annual", p$spatial.domain, ret, "rdata", sep=".") )
      if (file.exists( fn) ) load(fn)
      return ( BS )
    }

    if (exists( "libs", p)) RLibrary( p$libs )
    if ( is.null(ip)) ip = 1:p$nruns
     
    p0 = p # to make it explicit

    grids = unique( c(p0$spatial.domain.subareas , p0$spatial.domain ) ) # operate upon every domain
 
    for (gr in grids ) {
      print(gr)

      p1 = spatial_parameters( type=gr ) #target projection    
      L1 = bathymetry.db(p=p1, DS="baseline")

      tstatdir_p1 = file.path( project.datadirectory("bio.temperature"), "modelled", voi, p1$spatial.domain )

      O = array( NA, dim=c( nrow(L1), p0$ny, length(p$bstats)) )

      for ( iy in 1:p0$ny ) {

        y = p0$yrs[iy]

  			print ( paste("Year:", y)  )
        
  			V = temperature.db( p=p1, DS="predictions", yr=y, ret="sd"  )
  			V[ V < 0.1 ] = 100  # shrink weighting of unreasonably small SEs
  		  V[ which( !is.finite(V)) ] = 1000 # "
        W = 1/V^2   # weights: inverse variance, normalised
        W = W / rowSums(W)
        V = NULL
        
        P = temperature.db( p=p1, DS="predictions", yr=y, ret="mean"  )
        O[,iy,1] = c(apply( P*W, 1, sum, na.rm=T))
        O[,iy,2]  = c(apply( (P-rowMeans(P))^2*W, 1, sum, na.rm=T )) # weighted seasonal mean sums of squares
        W = NULL
  			O[,iy,3] = c(apply( P, 1, quantile, probs=p$lbm_quantile_bounds[1], na.rm=TRUE ))
        O[,iy,4] = c(apply( P, 1, quantile, probs=p$lbm_quantile_bounds[2], na.rm=TRUE ))
  			P = NULL
        O[,iy,5] = O[,iy,4] - O[,iy,3]
        O[,iy,6] = c(apply( P, 1, cumsum )) / p0$ny # normalized degree days (X 365 for degre days)
      }

      # save sp-time matrix for each stat .. easier to load into lbm this way   
      for ( st in 1:length(p$bstats) ){
        BS = O[,,st]
        fn = file.path( tstatdir_p1, paste("bottom.statistics.annual", p1$spatial.domain, p$bstats[st], "rdata", sep=".") )
        save( BS, file=fn, compress=T )
        BS = NULL
        gc()

      }

      # climatology
      clim = matrix( NA, nrow=nrow(L1), ncol=length(p$bstats) )
      for (si in 1:length(p$bstats)) {
        clim[,si] = rowMeans(O[,,si], na.rm=T)
      }
      fn.climatology = file.path( tstatdir_p1, paste("bottom.statistics.climatology", p1$spatial.domain, "rdata", sep=".") )
      colnames(clim) = p$bstats
      save( clim, file=fn.climatology, compress=T )
      clim = NULL
      gc()

    }

    return ("Completed")
  }


  # -----------------

  if (DS %in% c("spatial.annual.seasonal", "spatial.annual.seasonal.redo") ) {
    #\\ spatial, temporal (annual and seasonal) .. means only
    #\\ copy in array format for domain/resolution of interest for faster lookups
    outdir = file.path( project.datadirectory("bio.temperature"), "modelled", voi, p$spatial.domain )
    dir.create(outdir, recursive=T, showWarnings=F)

    outfile =  file.path( outdir, "temperature.spatial.annual.seasonal.rdata" )

    if ( DS=="spatial.annual.seasonal" ) {
      O = NULL
      if (file.exists(outfile)) load( outfile )
      return (O)
    }

    grids = unique( c( p$spatial.domain.subareas, p$spatial.domain) )
    for (gr in grids ) {
      #  print(gr)
      p1 = spatial_parameters( type=gr ) #target projection
      nlocs = nrow( bathymetry.db(p=p1, DS="baseline"))
      O = array( NA, dim=c(nlocs, p$ny, p$nw ) )
      for ( y in 1:p$ny ) {
        O[,y,] = temperature.db( p=p1, DS="predictions", yr=p$yrs[y], ret="mean" )
      }
      outdir = file.path( project.datadirectory("bio.temperature"), "modelled", voi, p1$spatial.domain )
      outfile =  file.path( outdir, "temperature.spatial.annual.seasonal.rdata" )
      save (O, file=outfile, compress=T )
      print(outfile)
    } 
    return( outfile )
  
  }


  # -----------------


  if (DS %in% c(  "timeslice", "timeslice.redo" )){

    tslicedir = file.path( project.datadirectory("bio.temperature"), "modelled", voi, p$spatial.domain )
    dir.create( tslicedir, showWarnings=F, recursive = TRUE )

    dyear_index = 1
    if (exists("dyears", p) & exists("prediction.dyear", p))  dyear_index = which.min( abs( p$prediction.dyear - p$dyears))

    if (DS %in% c("timeslice")) {
      O = NULL
      if (is.null(ret)) ret="mean"
      outfile =  file.path( tslicedir, paste("bottom.timeslice", dyear_index, ret, "rdata", sep=".") )
      if (file.exists( outfile ) ) load(outfile)
      return ( O )
    }

    p0 = p  # the originating parameters
    grids = unique( c( p$spatial.domain.subareas , p$spatial.domain) )

    for (gr in grids ) {
      # print(gr)
      p1 = spatial_parameters( type=gr ) #target projection
      tslicedir = file.path( project.datadirectory("bio.temperature"), "modelled", voi, p1$spatial.domain )
      nlocs = nrow( bathymetry.db(p=p1, DS="baseline"))
      O = matrix(NA, ncol=p$ny, nrow=nlocs)
      for ( r in 1:p$ny ) {
        P = temperature.db( p=p1, DS="predictions", yr=p$tyears[r], ret="mean"  )
        if (!is.null(P))  O[,r] = P[,dyear_index]
        P = NULL
      }
      outfileP =  file.path( tslicedir, paste("bottom.timeslice", dyear_index, "mean", "rdata", sep=".") )
      save( O, file=outfileP, compress=T )
      print(outfileP)

      O = NULL 
      
      O = matrix(NA, ncol=p$ny, nrow=nlocs)
      for ( r in 1:p$ny ) {
        V = temperature.db( p=p1, DS="predictions", yr=p$tyears[r], ret="sd"  )
        if (!is.null(V)) O[,r] = V[,dyear_index]
        V = NULL
      }
      outfileV =  file.path( tslicedir, paste("bottom.timeslice", dyear_index, "sd", "rdata", sep=".") )
      save( O, file=outfileV, compress=T )
      O = NULL
    }

    return ("Completed")
  }

  # ----------------------------


  if (DS %in% c(  "lbm.stats", "lbm.stats.redo" )){

    outdir = file.path(project.datadirectory("bio.temperature"), "modelled", voi, p$spatial.domain )
    
    if (DS %in% c("lbm.stats")) {
      stats = NULL
      fn = file.path( outdir, paste( "lbm.statistics", "rdata", sep=".") )
      if (file.exists(fn) ) load(fn) 
      return( stats )
    }

    # downscale and warp from p(0) -> p1
    # default domain
    S0 = lbm_db( p=p, DS="stats.to.prediction.grid" )
    Snames = colnames(S0)
    p0 = spatial_parameters( p=p, type=p$spatial.domain ) # from
    L0 = bathymetry.db( p=p0, DS="baseline" )
    L0i = lbm::array_map( "xy->2", L0[, c("plon", "plat")], gridparams=p0$gridparams )
    sreg = setdiff( p$spatial.domain.subareas, p$spatial.domain ) 

    for ( gr in sreg ) {
      p1 = spatial_parameters( p=p, type=gr ) # 'warping' from p -> p1
      L1 = bathymetry.db( p=p1, DS="baseline" )
      L1i = lbm::array_map( "xy->2", L1[, c("plon", "plat")], gridparams=p1$gridparams )
      L1 = planar2lonlat( L1, proj.type=p1$internal.crs )
      L1$plon_1 = L1$plon # store original coords
      L1$plat_1 = L1$plat
      L1 = lonlat2planar( L1, proj.type=p0$internal.crs )
      p1$wght = fields::setup.image.smooth( nrow=p1$nplons, ncol=p1$nplats, dx=p1$pres, dy=p1$pres, theta=p1$pres, xwidth=4*p1$pres, ywidth=4*p1$pres )
      stats = matrix( NA, ncol=ncol(S0), nrow=nrow(L1) )
      for ( i in 1:ncol(S0) ) {
        stats[,i] = spatial_warp( S0[,i], L0, L1, p0, p1, "fast", L0i, L1i )
      }
      colnames(stats) = Snames
      outdir_p1 = file.path(project.datadirectory("bio.temperature"), "modelled", voi, p1$spatial.domain)
      dir.create( outdir_p1, recursive=T, showWarnings=F )
      fn1_sg = file.path( outdir_p1, paste("lbm.statistics", "rdata", sep=".") )
      save( stats, file=fn1_sg, compress=T )
      print (fn1_sg)
    }
    return ("Completed")
  }



  # ----------------------------



  if (DS %in% c("complete", "complete.redo" )) {
    # static summaries
    
    if (DS=="complete") {
      TM = NULL
      outdir =  file.path( project.datadirectory("bio.temperature"), "modelled", voi, p$spatial.domain )
      outfile =  file.path( outdir, paste( "temperature", "complete", p$spatial.domain, "rdata", sep= ".") )
      if ( file.exists( outfile ) ) load( outfile )
      Tnames = names(TM)
      if (is.null(varnames)) varnames=Tnames
      varnames = intersect( Tnames, varnames )
      if (length(varnames) == 0) varnames=Tnames  # no match .. send all
      TM = TM[ , varnames]
      return(TM)
    }


    grids = unique( c(p$spatial.domain.subareas , p$spatial.domain ) ) # operate upon every domain
 
    for (gr in grids ) {
      print(gr)

      p1 = spatial_parameters( type=gr ) #target projection    
      L1 = bathymetry.db(p=p1, DS="baseline")

      BS = temperature.db( p=p1, DS="lbm.stats" )
      colnames(BS) = paste("t", colnames(BS), sep=".")
      TM = cbind( L1, BS )

      CL = temperature.db( p=p1, DS="bottom.statistics.climatology" )
      colnames(CL) = paste(p$bstats, "climatology", sep=".")
      TM = cbind( TM, CL )

      # bring in last stats
      outdir = file.path(project.datadirectory("bio.temperature"), "modelled", voi, p1$spatial.domain)
      dir.create( outdir, recursive=T, showWarnings=F )
      outfile =  file.path( outdir, paste( "temperature", "complete", p1$spatial.domain, "rdata", sep= ".") )
      save( TM, file=outfile, compress=T )

      print( outfile )

    }

    return( outdir )
  }

}
AtlanticR/bio.temperature documentation built on May 28, 2019, 11:35 a.m.