Removed prediction years and returnDataType

noaa-nwfsc · Nov 1, 2024 · 7d3f832 · 7d3f832
1 parent ac3ecc5
commit 7d3f832
Show file tree

Hide file tree

Showing 7 changed files with 15 additions and 54 deletions.
diff --git a/R/createStoredObjects.R b/R/createStoredObjects.R
@@ -36,14 +36,12 @@ min.lat=-90
 max.lat=90
 years=seq(1967, 2023, 1)
 months=seq(1,12,1)
-removeBering=FALSE
 
 # Function defined in create_OceanData_Object.R
 oceanData_ERSST<-getOceanData(dataSet=dataSet,
                         min.lon=min.lon, max.lon=max.lon,
                         min.lat=min.lat, max.lat=max.lat,
-                        years = years, months = months,
-                        removeBering = removeBering)
+                        years = years, months = months)
 save(x = "oceanData_ERSST", file = 'data/oceanSSTData.RData')
 load('data/oceanSSTData.RData')
 
@@ -70,13 +68,11 @@ min.lat=-90
 max.lat=90
 years=seq(1980, 2023, 1)
 months=seq(1,12,1)
-removeBering=FALSE
 
 oceanData_SSH<-getOceanData(dataSet=dataSet,
                             min.lon=min.lon, max.lon=max.lon,
                             min.lat=min.lat, max.lat=max.lat,
-                            years = years, months = months,
-                            removeBering = removeBering)
+                            years = years, months = months)
 save(x = "oceanData_SSH", file = 'data/oceanSSHData.RData')
 load('data/oceanSSHData.RData')
 

diff --git a/R/crossValidation.R b/R/crossValidation.R
@@ -1,5 +1,4 @@
 # Run LOO cross-validation, removing each year sequentially and predicting that year
-# This version uses the getIndex.R script to simplify the code
 
 #*************************************************************
 #  Create a loop for cross-validation
@@ -29,8 +28,7 @@ LOO_CV <- function(response = response,
                                min.lon = min.lon, max.lon = max.lon,
                                min.lat = min.lat, max.lat = max.lat,
                                years = years, years.fit = years.fit,
-                               months = months,
-                               removeBering = removeBering)
+                               months = months)
 
 
     # output

diff --git a/R/getOceanData.R b/R/getOceanData.R
@@ -21,15 +21,11 @@ library(ncdf4)
 # Create the function
 #***************************************************************
 getOceanData<-function(dataSet='ERSST',
-                 returnDataType='anom', returnObjectType='array',
                  min.lon=158, max.lon=246,
                  min.lat=10, max.lat=62,
-                 years=seq(1980, 2020, 1), months=seq(1,12,1),
-                 removeBering=TRUE) {
+                 years=seq(1980, 2020, 1), months=seq(1,12,1)) {
 
   # Check conditions
-  if (!returnDataType %in% c('anom','raw')) return (cat("returnDataType must be either 'anom' or 'raw'"))
-  if (!returnObjectType %in% c('array')) return (cat("returnObjectType must be 'array'"))
   if (min.lon < 0 | min.lon > 360 | max.lon < 0 | max.lon > 360) return (cat("Longitude is in degrees east (Tokyo is 139.8, Seattle is 237.6)"))
   if (min.lat < -90 | min.lat > 90 | max.lat < -90 | max.lat > 90) return (cat("Latitude is in degrees from the equator (should not exceed 90)"))
 
@@ -86,10 +82,7 @@ getOceanData<-function(dataSet='ERSST',
     dimnames(sst) <- list(seq(0,358,2), seq(-88,88,2), dates)
     # get anomaly, if requested
     for (mo in 1:nrow(year_mo)) {
-      if (returnDataType == 'anom') {
-          sst.temp <- sst[,, dates == year_mo$label[mo]] - sst.ltm[,, year_mo$month[mo]]
-      } else sst.temp <- sst[,, dates == year_mo$label[mo]]
-      returnData[, , mo]<-sst.temp
+      returnData[, , mo] <- sst[,, dates == year_mo$label[mo]] - sst.ltm[,, year_mo$month[mo]]
     }
     nc_close(ncin)
     returnData <- returnData[,89:1,] # I think this would fix the reversed latitude
@@ -103,27 +96,11 @@ getOceanData<-function(dataSet='ERSST',
         ssh.temp <- ncvar_get(ncin,"sshg")
       }
       # get anomaly, if requested
-      if (returnDataType == 'anom') {
-        ssh.temp[,, year_mo$month[ym]]<-ssh.temp[,, year_mo$month[ym]] - ssh.ltm[,, year_mo$month[ym]]
-      }
+      ssh.temp[,, year_mo$month[ym]]<-ssh.temp[,, year_mo$month[ym]] - ssh.ltm[,, year_mo$month[ym]]
       returnData[,,ym]<-ssh.temp[lon.index, lat.index, year_mo$month[ym]]
       if (year_mo$month[ym]==12) nc_close(ncin)
     }
   }
 
-  if (removeBering) {
-    # We extracted SST data for the full grid, but we don't want some portions of it
-    #  Remove the Bering Sea
-    returnData[lon.subset < 206, lat.subset > 56,] <- NA
-    returnData[lon.subset < 202, lat.subset > 54,] <- NA
-    returnData[lon.subset < 196, lat.subset > 52,] <- NA
-    returnData[lon.subset < 158, lat.subset > 50,] <- NA
-    returnData[lon.subset < 156, lat.subset > 48,] <- NA
-    returnData[lon.subset < 154, lat.subset > 46,] <- NA
-    returnData[lon.subset < 152, lat.subset > 44,] <- NA
-    returnData[lon.subset < 148, lat.subset > 42,] <- NA
-    returnData[lon.subset < 146, lat.subset > 40,] <- NA
-  }
-
   return(returnData)
 }
diff --git a/R/get_index.R b/R/get_index.R
@@ -7,8 +7,7 @@ get_CMISST_index <- function(response, oceanData=oceanData_ERSST,
                              years=NA, years.fit=years.fit,
                              months=1:12, years.pred=NA,
                              min.lon=158, max.lon=246,
-                             min.lat=10, max.lat=62,
-                             removeBering=TRUE) {
+                             min.lat=10, max.lat=62) {
   # Verify that the response is what we expect
   if (ncol(response)!=2) { print("incorrect data - requires a 2-column data frame with year and the response"); return(NA) }
   colnames(response)<-c("year","val")

diff --git a/R/makePlots.R b/R/makePlots.R
@@ -63,7 +63,7 @@ makeBiplot <- function(input.season = input.season, cmisst = cmisst) {
        x = par("usr")[1]*0.8, y=par("usr")[4]*0.80, cex=1.6, col="blue")
   if (input.loocv) {
     mae <- cmisst[[7]]
-    text(paste("MAE =", round(mae[season,"mae"], 2)),
+    text(paste("MAE =", round(mean(abs(mae[mae$season==input.season,"mae"])), 2)),
          x = par("usr")[1]*0.75, y=par("usr")[4]*0.60, cex=1.6, col="blue")
   }
 }
@@ -96,7 +96,8 @@ makeTimeSeriesPlot <- function(input.season = input.season, cmisst = cmisst,
   preds_new<-reverse_scale(preds_new, attr(response.tmp$val.scl, "scaled:center"), attr(response.tmp$val.scl, "scaled:scale"))
   if (input.log) preds_new<- exp(preds_new)
   # replace just the ones that were not used during fitting
-  preds[index$year %in% input.years.pred,]<-preds_new[index$year %in% input.years.pred,]
+  #preds[index$year %in% input.years.pred,]<-preds_new[index$year %in% input.years.pred,]
+  preds[is.na(index$val),]<-preds_new[is.na(index$val),]
 
   preds<-data.frame(preds)
   # unlag the year to show the plot in return year
@@ -136,7 +137,7 @@ makeLOOplot <- function(cmisst = cmisst, season = "spr") {
   index <- cmisst[[7]] # this is just the loo results
   index2<-index[index$season==season & index$model=="cmisst",]
   lines(index2$year, index2$pred, lwd=3, col="deepskyblue2")
-  text(labels = paste("LOO MAE CMISST =", round(mean(index2$mae),2)),
+  text(labels = paste("LOO MAE CMISST =", round(mean(abs(index2$mae)),2)),
        x = par("usr")[1]+9, y=par("usr")[4]*0.80, cex=1.0, col="deepskyblue2")
 }
 

diff --git a/R/run_CMISST.R b/R/run_CMISST.R
@@ -34,7 +34,6 @@ source('R/crossValidation.R')
 
 # set parameters (for most users, leave these as they are)
 months=seq(1,12,1)
-removeBering=FALSE
 
 # Parameters --------------------------------------------------------------
 
@@ -65,11 +64,6 @@ input.long= c(158, 246)
 #  Years after the most recent year in the response will be predicted
 input.years= c(1980, 2023)
 
-# Prediction years (ocean years)
-#  These years will not be included in calculating the CMISST index,
-#  but will be in the index output for use in a predictive model
-#input.years.pred=c(2020)
-input.years.pred=NA
 
 #************************************
 #  For Leave One Out Cross-validation

diff --git a/R/updateCMISST.R b/R/updateCMISST.R
@@ -14,8 +14,7 @@ updateCMISST <- function() {
   min.lat = input.lat[1]
   max.lat = input.lat[2]
   years = seq(input.years[1], input.years[2], 1)
-  years = sort(unique(c(years, input.years.pred)))
-
+
   response.tmp <- response
   # Lag the response variable
   response.tmp$year <- response.tmp$year - as.numeric(input.lag)
@@ -31,18 +30,15 @@ updateCMISST <- function() {
   response.tmp$val.scl <- scale(response.tmp$val)
 
   # Which years are being fit to?
-  if (!is.na(input.years.pred[1])) {
-    years.fit<-years[!years %in% input.years.pred & years %in% response.tmp$year] # will be needed to calculate the covariance
-  } else years.fit <- years[years %in% response.tmp$year]
+  years.fit <- years[years %in% response.tmp$year]
 
   # Calculate the CMISST index
   cmisst <- get_CMISST_index(response = response.tmp[,c("year","val.scl")],
-                             oceanData = oceanData, years.pred = input.years.pred,
+                             oceanData = oceanData, 
                              min.lon = min.lon, max.lon = max.lon,
                              min.lat = min.lat, max.lat = max.lat,
                              years = years, years.fit = years.fit,
-                             months = months,
-                             removeBering = removeBering)
+                             months = months)
 
   if (input.loocv) {
     loocv <- LOO_CV(response = response.tmp[,c("year","val.scl")],