############################################################################################ ## R Climate Script to analyze Mean Sea Level Anomaly Trends ## ## ## Download and process Netcdf file on mean sea level anomalies ## ## Source file: http://ibis.grdl.noaa.gov/SAT/SeaLevelRise/slr/slr_sla_pac_free_all_66.nc ## ## ## Developed by D Kelly O'Day 3/8/10 ## ## ## ## Data file is netcdf format. Use ndcf package to open and read data file ## ## Had to download file and save loacally to be able to use package(ncdf) ## ## Only altimetry measurements between 66°S and 66°N have been processed. ## ## An inverted barometer has been applied to the time series. ## ## The estimates of sea level rise do not include glacial isostatic adjustment effects ## ## on the geoid, which are modeled to be +0.2 to +0.5 mm/year when globally averaged. ## ############################################################################################ # Load libraries library(ncdf) par(ps=10); par(las=1); par(oma=c(2.5,2,1,1)); par(mar=c(1,4,2,1)) # Set local PC link, read netcdf (nc) file link <- "C:\\R_Home\\Charts & Graphs Blog\\RClimateTools\\NETCDF\\slr_sla_gbl_free_all_66.nc" nc1 <- open.ncdf(link) # print out netcdf file properties ## nc_data <- print(nc1) ## define series names, sizes, start-stop dates series <- c("sla_tx", "sla_e2", "sla_g1", "sla_j1", "sla_n1", "sla_j2") sat_name <- c("TOPEX", "ERS-2", "GFO", "Jason-1", "Envisat", "Jason-2") time<- c("time_tx", "time_e2", "time_g1", "time_j1", "time_n1", "time_j2") min_yr <- c(1992, 1995,2000, 2002, 2002,2008) min_mo <- c( 12, 5, 1, 1, 9, 7) #### Function to read and process each data series seperately ####### data_func <- function(x) { # Read Series data - get data from netcdf file msl <- get.var.ncdf(nc1, series[x]) t <- get.var.ncdf(nc1, time[x]) yr <- as.integer(t) mo <- ceiling((t-yr)/(0.083)) dt_c <- paste(yr,mo, 15, sep="/") dt <- as.Date(dt_c, format= "%Y/%m/%d") ### create data frame df_1 <- data.frame(t,yr, mo,dt, msl) ### Calc monthly avg msl mo_msl <- tapply(df_1$msl, INDEX = df_1$dt, FUN = mean, na.rm=T) mo_msl <- ts(mo_msl, start=c(min_yr[x], min_mo[x]), freq=12) dt <- time(mo_msl)+0.001 # adjust because Jan shows up as x.0 dt_yr <- as.integer(dt) dt_mo <- ceiling((dt -dt_yr)/(0.083)) df_ts <- data.frame(dt,dt_mo, dt_yr, mo_msl) names(df_ts) <- c("dt","mo", "yr", "msl") return(df_ts) } ###################################################################### # Process satellite series: # There are 6 series with time coverage overlap. Only need 3 series to get full time coverage # [1]-Topex, [4]- Jason-1,[6]- Jason-2. arrow_y <- c(-25,-24, -24, -22, -20, -20) # set y position for plot or sat coverage arrows ### Create Blank Plot note_1 <- "Altimetry data provided by the NOAA Laboratory for Satellite Altimetry\nGlobal ocean (60S to 60N)\nSeasonal Signal Removed, Inverted baromter applied." note_2 <- "Estimates do not include glacial isostatic adjustment (GIA) effects \nGIA estimated to be +0.2 to +0.5 mm/yr" msl_title <- "Mean Sea Level Trend: 1993 to 2010" plot(seq(1992, 2011), seq(1992,2011), type="n", col = "blue", xlim =c(1992,2011), ylim=c(-30,45),xlab="", ylab = "Mean Sea Level Anomaly - mm", main = msl_title, cex.lab = 0.85, cex.axis=0.85 ) text(1992, 40, note_1, cex=0.7, adj=0, col = "black") text(1992,34.7, note_2, cex=0.7,adj=0,col = "black") # Topex Data Series top_df <- data_func(1) top_df <- subset(top_df, top_df$dt < 2002.1) points(top_df$dt, top_df$msl, type="l", col = "black") t_min <- min(top_df$dt) t_max <- max(top_df$dt) t_avg <- mean(c(t_min, t_max)) arrows(t_min, arrow_y[1], t_max, arrow_y[1], length=0.025, code=3, col = "black") rect(t_avg-1.25, arrow_y[1]-0.5, t_avg+1.25, arrow_y[1]+0.5, col = "white", border = "white") text(t_avg, arrow_y[1], sat_name[1], cex = 0.8, col = "black") # Jason 1 Data Series j1_df <- data_func(4) j1_df <- subset(j1_df, j1_df$dt>2002.1 & j1_df$dt < 2008.58) points(j1_df$dt, j1_df$msl, type = "l", col = "blue") t_min <- min(j1_df$dt) t_max <- max(j1_df$dt) t_avg <- mean(c(t_min, t_max)) arrows(t_min, arrow_y[4], t_max, arrow_y[4], length=0.025, code=3, col = "blue") rect(t_avg-1.25, arrow_y[4]-0.5, t_avg+1.25, arrow_y[4]+0.5, col = "white", border = "white") text(t_avg,arrow_y[4], sat_name[4], cex = 0.8, col = "blue") # Jason 2 Data Series j2_df <- data_func(6) j2_df <- subset(j2_df, j2_df$dt >2008.58) points(j2_df$dt, j2_df$msl, type = "l", col = "green") t_min <- min(j2_df$dt) t_max <- max(j2_df$dt) t_avg <- mean(c(t_min, t_max)) arrows(t_min, arrow_y[6], t_max, arrow_y[6], length=0.025, code=3, col = "green") rect(t_avg-1.25,arrow_y[6]-0.5, t_avg+1.25,arrow_y[6] +0.5, col = "white", border = "white") text(t_avg,arrow_y[6], sat_name[6], cex = 0.8, col = "green") ## merge() 3 series - perform linear regression to get effective trend rate msl_df <- merge(top_df, j1_df,all=T) msl_df <- merge(msl_df, j2_df, all=T) last_row <- nrow(msl_df) last_mo <- msl_df[last_row, 2] last_yr <- msl_df[last_row,3] last_val <- msl_df[last_row, 4] last_date_c <- paste(last_mo, "/", 15, "/", last_yr, sep="") last_date_d <- as.Date(last_date_c, "%m/%d/%Y") last_month <- format(last_date_d, "%b") date_note <- paste(last_month, ", ", last_yr, " @ ", signif(last_val,3)," mm", sep="") points(msl_df[last_row,1], last_val, pch = 16, col="green") points(1993, 18, pch=16, col="green") text(1994, 18, date_note, cex=0.8, adj=0) ## Linear trend line msl_lm <- lm(msl_df$msl ~ msl_df$dt) lm_note <- paste("MSL rise - ", signif(coef(msl_lm)[2], 2), " mm per year") text(1994, 20, lm_note, adj = 0, cex = 0.85) a <- coef(msl_lm)[1] b <- coef(msl_lm)[2] ## Determine min & max yr_frac yr1 <- min(msl_df$dt) yr2 <- max(msl_df$dt) ## calc y vals y1 <- a + b* yr1 y2 <- a + b *yr2 ## set x_val & y_val for start & end points of rgression line x_val <- c(yr1, yr2) y_val <- c(y1, y2) # add lines to plot lines(x_val,y_val, type = "l", col = "red") b <- signif(b, 3) points(x=c(1992,1993.5), y=c(20,20),col="red", type = "l") ## Plot Annotation my_date <- format(Sys.time(), "%m/%d/%y") mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", side = 1, line = 1, cex=0.7, outer = T, adj = 0) mtext(my_date, side = 1, line =1, cex = 0.7, outer = T, adj = 1)