
# MAKE DATASETS OF PALEO-CLIMATE and NPP
# from the cleaned dataset of paleo-climatic records (Data_Paleoclimate_Regions.txt)

# continental (NA, EU) and regional (6 regions)
# with global time bins (stratigraphic stages)

# Contents
# line  30: cunningly sum paleoclimatic localities into grid cells
# line 375: calculate NPP for each grid cell
# line 450: make regional and continental datasets from the gridded ones
#           and add the information into Data_Regions and Data_Continents

# produces: 
# Data_GriddedPaleoclimate (intermediate file, not included in upload)

# modifies Data_Regions and Data_Continents as included in upload
# (adds NPP data to the mammal data generated in Rscript_MammalDataProcessing.r)


# Script by Susanne Fritz, last changed 6th Oct 2015

#############################################################################


# sum the paleoclimatic information into one-degree cells
# and calculate NPP from the climate dataset

# function to classify one min and max value into a range of bins
# minval, maxval - numeric values, the min and max value between which you want to classify
# rangedat - a vector of break points for the bins to classify
# returns a vector of TRUE/FALSE matching rangedat,
# which tells you which bins the data between min and max values fall into
makebintable <- function(minval, maxval, rangedat) {
  ret <- data.frame(Bins=rangedat, Min=minval<=rangedat, Max=maxval>=rangedat)
  ret$both <- ret$Min & ret$Max
  return(as.numeric(ret$both))
}

# wrapper function to use makebintable on several min and max values
# and get out a frequency distribution
# minx, maxx - vectors of min and max values
# steps - specification of the bin size to use
# returns a vector of the bin values and
# how often they were counted in between the min and max values
binmydata <- function(minx, maxx, steps=0.1) {
  
  # first, remove NAs
  isna <- which(is.na(minx) | is.na(maxx))
  if (length(isna)>0) {
    minx <- minx[-isna]
    maxx <- maxx[-isna]
  }
  # then, get the extremes and make the bins using the "steps" argument
  extremes <- c(floor(min(minx)), ceiling(max(maxx)))
  rangedat <- seq(extremes[1], extremes[2], by=steps)
  # get a vector or matrix that says which bin each data point or range falls into
  for (i in seq(along=minx)) {
    linei <- makebintable(minx[i], maxx[i], rangedat)
    if (i==1) (outp <- linei) else (outp <- rbind(outp, linei))
  }
  # if there are several data points, count the frequency for each bin
  myeach <- colSums(outp)
  rangedat <- rangedat[myeach!=0]
  myeach <- myeach[myeach!=0]
  # return a vector with bin values repeated as often as they were counted
  return(rep(rangedat, times=myeach))
}


setwd("D://Data//Gammadiv//FilesForUpload")
library(hdrcde)
# load the paleoclimatic dataset
cl <- read.delim("Data_Paleoclimate_Regions.txt", as.is=TRUE)
# this has a row for each location-stage combination
head(cl)
# data checking
all(is.na(cl$min_MAT_C[is.na(cl$max_MAT_C)]))
all(is.na(cl$max_MAT_C[is.na(cl$min_MAT_C)]))
all(is.na(cl$min_MAP_mm_yr[is.na(cl$max_MAP_mm_yr)]))
all(is.na(cl$max_MAP_mm_yr[is.na(cl$min_MAP_mm_yr)]))
cl$isna <- is.na(cl$min_MAT_C) & is.na(cl$min_MAP_mm_yr)
summary(as.factor(cl$Source))
summary(cl$WM_ID)
summary(as.factor(cl$AutomaticStage)) # 423 lines for automatic dates
summary(as.factor(cl$BestStage)) # 439 lines for best dates
luni <- function(x) (length(unique(x)))
tapply(cl$Data_ID[which(is.na(cl$BestStage))], cl$Source[which(is.na(cl$BestStage))], luni)
tapply(cl$Spatial_ID[which(is.na(cl$BestStage))], cl$Source[which(is.na(cl$BestStage))], luni)
tapply(cl$Data_ID[which(!is.na(cl$BestStage))], cl$Source[which(!is.na(cl$BestStage))], luni)
tapply(cl$Spatial_ID[which(!is.na(cl$BestStage))], cl$Source[which(!is.na(cl$BestStage))], luni)
tapply(cl$MAT_C, cl$WM_ID, length) # some grid cells have a lot of points in them

# 50 entries without min and max, just a point estimate for MAT
length(which(!is.na(cl$MAT_C) & (is.na(cl$MAT_error) | cl$MAT_error==0))) 
# check the range for points where min and max are given
cl$MAT_range <- cl$max_MAT_C - cl$min_MAT_C
summary(cl$min_MAT_C)
summary(cl$MAT_range[!is.na(cl$MAT_error) & cl$MAT_error>0])
cl[which(cl$MAT_error>10),]
# one data point had the same error as mean - probably a mistake!
# set the error to NA and entered the estimate as min and max as well
cl[which(cl$MAT_C==cl$MAT_error),]
cl$min_MAT_C[cl$Data_ID==25] <- cl$MAT_C[cl$Data_ID==25]
cl$max_MAT_C[cl$Data_ID==25] <- cl$MAT_C[cl$Data_ID==25]
cl$MAT_error[cl$Data_ID==25] <- NA
cl[which(cl$Data_ID==25),]
summary(cl$MAT_range[!is.na(cl$MAT_error) & cl$MAT_error>0])
# for MAP: 38 entries without min and max
length(which(!is.na(cl$MAP_mm_yr) & (is.na(cl$MAP_error) | cl$MAP_error==0))) 
# check the range for points where min and max are given
cl$MAP_range <- cl$max_MAP_mm_yr - cl$min_MAP_mm_yr
summary(cl$min_MAP_mm_yr)
summary(cl$MAP_range[!is.na(cl$MAP_error) & cl$MAP_error>0])
cl[which(cl$MAP_error>500),] # these look fine
summary(cl) # check number of NAs are identical for the mean, min and max columns (error can have more NAs, this should be covered)
# get numbers (not included in paper text as NPP calculation in the grid
# cells further reduced the numbers):
summary(as.factor(cl$Continent))
summary(as.factor(cl$Region))
length(unique(cl$Spatial_ID))
tapply(cl$Spatial_ID, cl$Continent, luni)
summary(as.factor(cl$BestStage))
# how many locations have both MAT and MAP? 
cl$bothna <- is.na(cl$MAT_range) | is.na(cl$MAP_range)
summary(cl$bothna)
length(cl[which(!cl$bothna & is.na(cl$BestStage)),1]) # 357 entries for automatic dating...
length(cl[which(!cl$bothna & !is.na(cl$BestStage)),1]) # 373 entries for best dates
length(unique(cl$Data_ID[!cl$bothna & is.na(cl$BestStage)])) # 228 records for automatic dating
length(unique(cl$Data_ID[!cl$bothna & !is.na(cl$BestStage)])) # 292 records for best dating
length(unique(cl$Spatial_ID[!cl$bothna & is.na(cl$BestStage)])) # 109 locations for automatic dating
length(unique(cl$Spatial_ID[!cl$bothna & !is.na(cl$BestStage)])) # 136 locations for automatic dating
# entries per continent, automatic dates
tapply(cl$Data_ID[!cl$bothna & is.na(cl$BestStage)], cl$Continent[!cl$bothna & is.na(cl$BestStage)], length)
# entries per continent, best dates
tapply(cl$Data_ID[!cl$bothna & !is.na(cl$BestStage)], cl$Continent[!cl$bothna & !is.na(cl$BestStage)], length)
# locations per continent, automatic dates
tapply(cl$Spatial_ID[!cl$bothna & is.na(cl$BestStage)], cl$Continent[!cl$bothna & is.na(cl$BestStage)], luni)
# locations per continent, best dates
tapply(cl$Spatial_ID[!cl$bothna & !is.na(cl$BestStage)], cl$Continent[!cl$bothna & !is.na(cl$BestStage)], luni)
# grid cells per continent, automatic dates
tapply(cl$WM_ID[!cl$bothna & is.na(cl$BestStage)], cl$Continent[!cl$bothna & is.na(cl$BestStage)], luni)
# grid cells per continent, best dates
tapply(cl$WM_ID[!cl$bothna & !is.na(cl$BestStage)], cl$Continent[!cl$bothna & !is.na(cl$BestStage)], luni)


# loop through the grid cells and make a plant dataset for each
myids <- unique(cl$WM_ID) # 100 grid cells
length(myids) # 100 grid cells contain climatic data
# need to go through the possible data configurations:
# automatic stages and best stages
mydata <- c("AutoPlants", "BestPlants")

rm(output, newlines) # in case I repeat the loop!
# for visual checking, histograms for each grid cell are produced showing the climate data
pdf("Fig_HistogramsClimatePerCell.pdf", width=4, height=8)
par(mfrow=c(2,1))
for (i in seq(along=myids)) {
  print(i)
  # # which time bins are covered in any data in grid cell i?
  binsi <- unique(c(cl$AutomaticStage[which(cl$WM_ID==myids[i])], cl$BestStage[which(cl$WM_ID==myids[i])]))
  binsi <- binsi[!is.na(binsi)]
  for (j in seq(along=binsi)) { # loop through these bins
    # dataset for bin j in cell i
    cli <- cl[which(cl$WM_ID==myids[i] & (cl$AutomaticStage==binsi[j] | cl$BestStage==binsi[j])),]
    # loop through the 2 possible datasets
    for (k in seq(along=mydata)) {
      if (k==1) {
        cli1 <- cli[which(!is.na(cli$AutomaticStage)),]
      } else {
        cli1 <- cli[which(!is.na(cli$BestStage)),]
      }
      # only continue if the dataset contains rows
      if (length(cli1[,1])>0) {
      
        # data processing for grid cell i, time bin j, dataset k
        # first, MAT: see if we have any MAT data points
        if (any(!is.na(cli1$min_MAT_C))) {
          # if we only have one MAT data point, or if there are several data points that are too similar for the hdr function to work
          if (length(which(!is.na(cli1$min_MAT_C)))==1 | all(length(which(!is.na(cli1$min_MAT_C)))>1, length(unique(round(cli1$min_MAT_C[which(!is.na(cli1$min_MAT_C))],0)))==1 & length(unique(round(cli1$max_MAT_C[which(!is.na(cli1$min_MAT_C))],0)))==1)) {
            # make a plot for the one data point, showing the range or the point estimate like the frequency distribution for hdr
            if (length(which(!is.na(cli1$min_MAT_C)))==1) {
          	  plot(c(rep(cli1$min_MAT_C[!is.na(cli1$min_MAT_C)], 2), rep(cli1$max_MAT_C[!is.na(cli1$max_MAT_C)], 2)), c(0,1,1,0), type="l", xlab="tmp", ylab="Frequency", main=paste(mydata[k], " ID ", myids[i], ", ", binsi[j], ", 1 data point", sep=""))
            } else {
            	# if there are several very similar data points, plot the frequency distribution anyway as below
            	# i.e. make a vector of binned temperature values (fills bins of size 0.01 degrees using each point or range estimate)
            	tmp <- binmydata(cli1$min_MAT_C, cli1$max_MAT_C, steps=0.01)
          	  # then plot these
          	  hist(tmp, main=paste(mydata[k], " ID ", myids[i], ", ", binsi[j], ", ", length(which(!is.na(cli1$min_MAT_C))), " data points", sep=""))
            }
            # make up the mean from just the one min and max value
            meantmp <- mean(c(cli1$min_MAT_C, cli1$max_MAT_C), na.rm=TRUE)
            # because hdr does not work well with only one data point (ie one min and max),
            # needs some more cases here
            if (all(length(unique(round(cli1$min_MAT_C, 1)))==1 & length(unique(round(cli1$max_MAT_C, 1)))==1, unique(round(cli1$min_MAT_C,1))==unique(round(cli1$max_MAT_C,1)))) {
			         # if there is only one data point, and that is a point estimate,
      	       # set the standard deviation to 1 degree to get a reasonably high variance
	             # for weighting later
            	sdtmp <- 1
            	# make up the hdr object with the mean of min and max as the mode,
            	# and the appropriate values for 95th centile (2 standard deviations from the mean)
            	# and for the 50th centile (1 standard deviation from the mean)
            	hdri <- list(hdr=matrix(c(meantmp - sdtmp*3, meantmp + sdtmp*3, meantmp - sdtmp*2, meantmp + sdtmp*2, meantmp - sdtmp, meantmp + sdtmp), nrow=3, byrow=TRUE), mode=meantmp)
              # set the mean minimum value to the observed minus 2 standard deviations
              # to avoid minimum values that are higher than the mean and maxima that are lower
              mintmp <- meantmp - sdtmp*2
              maxtmp <- meantmp + sdtmp*2
            } else {
              # if there is one data point that is a range, or if there are several very similar data points,
            	# decide that the given range between min and max gives 4 units of standard deviation
            	sdtmp <- (max(cli1$max_MAT_C, na.rm=TRUE) - min(cli1$min_MAT_C, na.rm=TRUE))/4
              # make up the hdr object with the mean of min and max as the mode,
            	# and the appropriate values for 95th centile (the original min and max given)
            	# and for the 50th centile (1 standard deviation from the mean)
            	hdri <- list(hdr=matrix(c(meantmp - sdtmp*3, meantmp + sdtmp*3, min(cli1$min_MAT_C, na.rm=TRUE), max(cli1$max_MAT_C, na.rm=TRUE), meantmp - sdtmp, meantmp + sdtmp), nrow=3, byrow=TRUE), mode=meantmp)
              # in this case, min and max are the min and max of the original min and max values
              mintmp <- mean(cli1$min_MAT_C, na.rm=TRUE)
            	maxtmp <- mean(cli1$max_MAT_C, na.rm=TRUE)
            }
            vartmp <- sdtmp^2 # variance
            # make up the number of binned temperature values as 2 because I only used min and max
            sizetmp <- 2
          } else {
            # if we have more than one MAT data point, and they are quite different,
            # use the hdr function to get some more elaborate statistics:
            tmp <- binmydata(cli1$min_MAT_C, cli1$max_MAT_C, steps=0.01)
            hist(tmp, main=paste(mydata[k], " ID ", myids[i], ", ", binsi[j], ", ", length(which(!is.na(cli1$min_MAT_C))), " data points", sep=""))
            hdri <- hdr(tmp)
            # get mean, variance, and number of the binned MAT values
            meantmp <- mean(tmp)
            vartmp <- var(tmp)
            sizetmp <- length(tmp)
            mintmp <- mean(cli1$min_MAT_C, na.rm=TRUE)
            maxtmp <- mean(cli1$max_MAT_C, na.rm=TRUE)
          }
        } else (plot.new()) # if there are no MAT points, make an empty plot so the MAP plot can be underneath
        
        # same for precipitation: if we have any points...
        if (any(!is.na(cli1$min_MAP_mm_yr))) {
          # if we have one data point, or effectively one data pint
          if (length(which(!is.na(cli1$min_MAP_mm_yr)))==1 | all(length(which(!is.na(cli1$min_MAP_mm_yr)))>1, length(unique(round(cli1$min_MAP_mm_yr[which(!is.na(cli1$min_MAP_mm_yr))]),0))==1 & length(unique(round(cli1$max_MAP_mm_yr[which(!is.na(cli1$min_MAP_mm_yr))],0)))==1)) {
    		    # plotting
            if (length(which(!is.na(cli1$min_MAP_mm_yr)))==1) {
            	plot(c(rep(cli1$min_MAP_mm_yr[!is.na(cli1$min_MAP_mm_yr)], 2), rep(cli1$max_MAP_mm_yr[!is.na(cli1$max_MAP_mm_yr)], 2)), c(0,1,1,0), type="l", xlab="prc", ylab="Frequency", main=paste(mydata[k], " ID ", myids[i], ", ", binsi[j], ", 1 data point", sep=""))
            } else {
            	prc <- binmydata(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr, steps=1)
            	hist(prc, main=paste(mydata[k], " ID ", myids[i], ", ", binsi[j], ", ", length(which(!is.na(cli1$min_MAP_mm_yr))), " data points", sep=""))
            }
            # getting the values
            meanprc <- mean(c(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr), na.rm=TRUE)
            if (all(length(unique(round(cli1$min_MAP_mm_yr, 1)))==1 & length(unique(round(cli1$max_MAP_mm_yr, 1)))==1, unique(round(cli1$min_MAP_mm_yr,1))==unique(round(cli1$max_MAP_mm_yr,1)))) {
            	sdprc <- 120
            	hdrj <- list(hdr=matrix(c(meanprc - sdprc*3, meanprc + sdprc*3, meanprc - sdprc*2, meanprc + sdprc*2, meanprc - sdprc, meanprc + sdprc), nrow=3, byrow=TRUE), mode=meanprc)
            	minprc <- meanprc - sdprc*2
            	maxprc <- meanprc + sdprc*2
            } else {
            	sdprc <- (max(cli1$max_MAP_mm_yr, na.rm=TRUE) - min(cli1$min_MAP_mm_yr, na.rm=TRUE))/4
            	hdrj <- list(hdr=matrix(c(meanprc - sdprc*3, meanprc + sdprc*3, min(cli1$min_MAP_mm_yr, na.rm=TRUE), max(cli1$max_MAP_mm_yr, na.rm=TRUE), meanprc - sdprc, meanprc + sdprc), nrow=3, byrow=TRUE), mode=meanprc)
            	minprc <- mean(cli1$min_MAP_mm_yr, na.rm=TRUE)
            	maxprc <- mean(cli1$max_MAP_mm_yr, na.rm=TRUE)
            }
            varprc <- sdprc^2 # variance
            sizeprc <- 2
          # if we have several data points  
          } else {
		        prc <- binmydata(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr, steps=1)
      	    hist(prc, main=paste(mydata[k], " ID ", myids[i], ", ", binsi[j], ", ", length(which(!is.na(cli1$min_MAP_mm_yr))), " data points", sep=""))
            hdrj <- hdr(prc)
            meanprc <- mean(prc)
            varprc <- var(prc)
            sizeprc <- length(prc)
      	    minprc <- mean(cli1$min_MAP_mm_yr, na.rm=TRUE)
		        maxprc <- mean(cli1$max_MAP_mm_yr, na.rm=TRUE)
          }
        } else (plot.new()) # if there are no MAP points, make an empty plot so the next MAT plot can be on a new page
          
        # put it all together
        if (all(any(!is.na(cli1$min_MAT_C)), any(!is.na(cli1$min_MAP_mm_yr)))) {
          # if we have both MAT and MAP, make a line of data for grid cell i in time bin j for data type k
          newline <- data.frame(WM_ID=myids[i], Stage=binsi[j], MAT_NumData=length(which(!is.na(cli1$min_MAT_C))), MAT_hdr=mean(hdri$mode), MAT_95low=hdri$hdr[2,1], MAT_95high=max(hdri$hdr[2,], na.rm=TRUE), MAT_50low=hdri$hdr[3,1], MAT_50high=max(hdri$hdr[3,], na.rm=TRUE), MAT_meanofdist=meantmp, MAT_varofdist=vartmp, MAT_sizeofdist=sizetmp, MAT_meanofpoints=mean(c(cli1$min_MAT_C, cli1$max_MAT_C), na.rm=TRUE), MAT_varofpoints=var(c(cli1$min_MAT_C, cli1$max_MAT_C), na.rm=TRUE), MAT_meanmin=mintmp, MAT_meanmax=maxtmp, MAP_NumData=length(which(!is.na(cli1$min_MAP_mm_yr))), MAP_hdr=mean(hdrj$mode), MAP_95low=hdrj$hdr[2,1], MAP_95high=max(hdrj$hdr[2,], na.rm=TRUE), MAP_50low=hdrj$hdr[3,1], MAP_50high=max(hdrj$hdr[3,], na.rm=TRUE), MAP_meanofdist=meanprc, MAP_varofdist=varprc, MAP_sizeofdist=sizeprc, MAP_meanofpoints=mean(c(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr), na.rm=TRUE), MAP_varofpoints=var(c(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr), na.rm=TRUE), MAP_meanmin=minprc, MAP_meanmax=maxprc)
        } else { # if we only have MAT or MAP, make sure the other columns are set to NA
          if (any(!is.na(cli1$min_MAP_mm_yr))) { # only MAP data
            newline <- data.frame(WM_ID=myids[i], Stage=binsi[j], MAT_NumData=length(which(!is.na(cli1$min_MAT_C))), MAT_hdr=NA, MAT_95low=NA, MAT_95high=NA, MAT_50low=NA, MAT_50high=NA, MAT_meanofdist=NA, MAT_varofdist=NA, MAT_sizeofdist=NA, MAT_meanofpoints=NA, MAT_varofpoints=NA, MAT_meanmin=NA, MAT_meanmax=NA, MAP_NumData=length(which(!is.na(cli1$min_MAP_mm_yr))), MAP_hdr=mean(hdrj$mode), MAP_95low=hdrj$hdr[2,1], MAP_95high=max(hdrj$hdr[2,], na.rm=TRUE), MAP_50low=hdrj$hdr[3,1], MAP_50high=max(hdrj$hdr[3,], na.rm=TRUE), MAP_meanofdist=meanprc, MAP_varofdist=varprc, MAP_sizeofdist=sizeprc, MAP_meanofpoints=mean(c(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr), na.rm=TRUE), MAP_varofpoints=var(c(cli1$min_MAP_mm_yr, cli1$max_MAP_mm_yr), na.rm=TRUE), MAP_meanmin=minprc, MAP_meanmax=maxprc)
          } else { # only MAT data
            newline <- data.frame(WM_ID=myids[i], Stage=binsi[j], MAT_NumData=length(which(!is.na(cli1$min_MAT_C))), MAT_hdr=mean(hdri$mode), MAT_95low=hdri$hdr[2,1], MAT_95high=max(hdri$hdr[2,], na.rm=TRUE), MAT_50low=hdri$hdr[3,1], MAT_50high=max(hdri$hdr[3,], na.rm=TRUE), MAT_meanofdist=meantmp, MAT_varofdist=vartmp, MAT_sizeofdist=sizetmp, MAT_meanofpoints=mean(c(cli1$min_MAT_C, cli1$max_MAT_C), na.rm=TRUE), MAT_varofpoints=var(c(cli1$min_MAT_C, cli1$max_MAT_C), na.rm=TRUE), MAT_meanmin=mintmp, MAT_meanmax=maxtmp, MAP_NumData=length(which(!is.na(cli1$min_MAP_mm_yr))), MAP_hdr=NA, MAP_95low=NA, MAP_95high=NA, MAP_50low=NA, MAP_50high=NA, MAP_meanofdist=NA, MAP_varofdist=NA, MAP_sizeofdist=NA, MAP_meanofpoints=NA, MAP_varofpoints=NA, MAP_meanmin=NA, MAP_meanmax=NA)
          }
        }
        # transfer the information which data source was used into the right columns
        if (k==1) {
          names(newline)[2] <- "AutomaticStage"
          newline$BestStage <- NA
          newline <- newline[,c(1,2,29,3:28)]
        } else {
          names(newline)[2] <- "BestStage"
          newline$AutomaticStage <- NA
          newline <- newline[,c(1,29,2:28)]
        }
        # merge the newline datasets from different data sources together
        if (exists("newlines")) (newlines <- rbind(newlines, newline)) else (newlines <- newline)
      } # if the dataset type k for stage j and cell i contains rows
    } # loop through the possible datasets
    if (exists("output")) (output <- rbind(output, newlines)) else (output <- newlines)
    if (exists("newlines")) (rm(newlines))
  } # loop through the bins
} # loop through the cells
dev.off()
head(output)
any(duplicated(output[,1:3])) # there should be no duplicates for this
length(unique(output$WM_ID)) # should be 100
summary(output) # check no values are negative or something weird
summary(as.factor(output$AutomaticStage))
summary(as.factor(output$BestStage))
summary(as.factor(output$WM_ID))
# how many data points / grid cells did the different data sets end up with?
lengthnotna <- function(x) (length(x[!is.na(x)]))
luni <- function(x) (length(unique(x[!is.na(x)])))
# how many grid cells with just one data point?
tapply(output$AutomaticStage[which(output$MAT_NumData==1)], output$AutomaticStage[which(output$MAT_NumData==1)], lengthnotna)
tapply(output$AutomaticStage[which(output$MAP_NumData==1)], output$AutomaticStage[which(output$MAP_NumData==1)], lengthnotna)
tapply(output$BestStage[which(output$MAT_NumData==1)], output$BestStage[which(output$MAT_NumData==1)], lengthnotna)
tapply(output$BestStage[which(output$MAP_NumData==1)], output$BestStage[which(output$MAP_NumData==1)], lengthnotna)

# set varofpoints to 1 degree if it is 0 (no min and max given, just point estimate)
# because this is the variance for a standard deviation of 1
output[which(output$MAT_varofpoints==0),]
output$MAT_varofpoints[which(output$MAT_varofpoints==0)] <- 1
# how many cells have only this one point estimate?
summary(as.factor(output$BestStage[which(output$MAT_varofpoints==1 & output$MAT_NumData==1)]))
summary(as.factor(output$AutomaticStage[which(output$MAT_varofpoints==1 & output$MAT_NumData==1)]))
length(unique(output$WM_ID[output$MAT_varofpoints==1 & output$MAT_NumData==1]))
# check there are no weird things going on if there are
# no or just 1 data points in a cell
length(which(output$MAT_NumData==1))
length(unique(output$WM_ID[output$MAT_NumData==1]))
length(which(output$MAT_NumData==1 & !is.na(output$AutomaticStage)))
length(which(output$MAT_NumData==1 & !is.na(output$BestStage)))
length(unique(output$WM_ID[which(output$MAT_NumData==1 & !is.na(output$AutomaticStage))]))
length(unique(output$WM_ID[which(output$MAT_NumData==1 & !is.na(output$BestStage))]))
tapply(output$WM_ID[output$MAT_NumData==1], output$AutomaticStage[output$MAT_NumData==1], luni)
summary(output[which(output$MAT_NumData==0),])
summary(output[which(output$MAT_NumData==1),])
# for MAP, set it to 14400 (the variance for sd=120)
output[which(output$MAP_varofpoints==0),]
output$MAP_varofpoints[which(output$MAP_varofpoints==0)] <- 14400
summary(as.factor(output$BestStage[which(output$MAP_varofpoints==14400 & output$MAP_NumData==1)]))
summary(as.factor(output$AutomaticStage[which(output$MAP_varofpoints==14400 & output$MAP_NumData==1)]))
length(unique(output$WM_ID[output$MAP_varofpoints==14400 & output$MAP_NumData==1]))
# check there are no weird things going on if there are
# no or just 1 data points in a cell
length(which(output$MAP_NumData==1))
length(which(output$MAP_NumData==1 & !is.na(output$AutomaticStage)))
length(which(output$MAP_NumData==1 & !is.na(output$BestStage)))
length(unique(output$WM_ID[output$MAP_NumData==1]))
length(unique(output$WM_ID[which(output$MAP_NumData==1 & !is.na(output$AutomaticStage))]))
length(unique(output$WM_ID[which(output$MAP_NumData==1 & !is.na(output$BestStage))]))
tapply(output$WM_ID[output$MAP_NumData==1], output$AutomaticStage[output$MAP_NumData==1], luni)
tapply(output$WM_ID[output$MAP_NumData==1], output$BestStage[output$MAP_NumData==1], luni)
summary(output[which(output$MAP_NumData==0),])
summary(output[which(output$MAP_NumData==1),])

# append missing information
head(cl)
output$Region <- cl$Region[match(output$WM_ID, cl$WM_ID)]
output$Continent <- cl$Continent[match(output$WM_ID, cl$WM_ID)]
table(output$Region, output$BestStage)
table(output$Region, output$AutomaticStage)
table(output$Continent, output$BestStage)
table(output$Continent, output$AutomaticStage)
write.table(output, "Data_GriddedPaleoclimate.txt", sep="\t", row.names=FALSE, col.names=TRUE)


#############################################################################

# calculate NPP following the Miami model for the grid cells with both
# MAT and MAP data
# further modifies the Data_GriddedPaleoclimate.txt file generated above.

setwd("D://Data//Gammadiv//FilesForUpload")
lengthnotna <- function(x) (length(x[which(!is.na(x))])) # function to count the number of non-NA entries in a vector
luni <- function(x) (length(unique(x))) # function to count unique entries

# calculate NPP as per the Miami model
# - NPP based on precipitation if this NPP value is smaller than
# the NPP value based on temperature; else NPP based on temperature
# - returns NA if one of the two input values is NA
npp <- function(mat, map) {
	nppT <- 3000/(1+exp(1.315-0.119*mat))
	nppP <- 3000*(1-exp(-0.000664*map))
	NPP <- nppT
	NPP[which(nppP < nppT)] <- nppP[which(nppP < nppT)]
	NPP[which(is.na(nppT) | is.na(nppP))] <- NA
	return(NPP)
}
# get which value the calculation was based on
nppwhich <- function(mat, map) {
  nppT <- 3000/(1+exp(1.315-0.119*mat))
  nppP <- 3000*(1-exp(-0.000664*map))
  whichnpp <- rep("nppT", length(nppT))
  whichnpp[which(nppP < nppT)] <- "nppP"
  whichnpp[which(is.na(nppT) | is.na(nppP))] <- NA
  return(whichnpp)
}
# calculate NPP as per the Miami model BUT user selects which
# of the two NPP values is used - the one based on temperature or
# the one based on precipitation
# - needed this for the confidence intervals, which are otherwise 
# sometimes based on temperature for a mean that is based on precip
# and vice versa
npp.select <- function(mat, map, whichnpp) {
  nppT <- 3000/(1+exp(1.315-0.119*mat))
  nppP <- 3000*(1-exp(-0.000664*map))
  NPP <- nppT
  NPP[which(whichnpp=="nppP")] <- nppP[which(whichnpp=="nppP")]
  NPP[which(is.na(nppT) | is.na(nppP) | is.na(whichnpp))] <- NA
  return(NPP)
}

# run the calculations to get NPP
cli <- read.delim("Data_GriddedPaleoclimate.txt", as.is=TRUE)
head(cli)
cli <- cli[,c(1,30,31,2:29)]
names(cli)
cli$NPP_hdr <- npp(cli$MAT_hdr, cli$MAP_hdr) # the highest density point (mean estimate)
cli$NPP_which <- nppwhich(cli$MAT_hdr, cli$MAP_hdr) # which climate variable was final NPP derived from?
cli$NPP_50low <- npp.select(cli$MAT_50low, cli$MAP_50low, cli$NPP_which) # the lower bound of the 50% interval (for weighting by the interval breadth later)
cli$NPP_50high <- npp.select(cli$MAT_50high, cli$MAP_50high, cli$NPP_which)
cli$NPP_meanmin <- npp.select(cli$MAT_meanmin, cli$MAP_meanmin, cli$NPP_which) # the mean of the minimum estimates
cli$NPP_meanmax <- npp.select(cli$MAT_meanmax, cli$MAP_meanmax, cli$NPP_which)
# investigate the result
summary(cli$NPP_hdr) # calculated NPP for 290 cell-stage combinations in total
tapply(cli$NPP_hdr, cli$AutomaticStage, summary)
length(which(!is.na(cli$AutomaticStage))) # 186 cell-stage combinations
length(which(!is.na(cli$BestStage))) # 194 cell-stage combinations
# how many grid cells for each dataset that have an NPP estimate?
luni(cli$WM_ID[which(is.na(cli$AutomaticStage) & !is.na(cli$NPP_hdr))])
luni(cli$WM_ID[which(is.na(cli$BestStage) & !is.na(cli$NPP_hdr))])
# how many cell-stage combinations that have an NPP estimate?
lengthnotna(cli$NPP_hdr[which(!is.na(cli$AutomaticStage))]) # 140 combinations with NPP data
lengthnotna(cli$NPP_hdr[which(!is.na(cli$BestStage))]) # 150
summary(cli$NPP_50low)
summary(cli$NPP_50high)
summary(cli$NPP_meanmin)
summary(cli$NPP_meanmax)
cli[which(cli$Continent=="North America" & cli$BestStage=="Early Pliocene"),]
head(cli)
write.table(cli, "Data_GriddedPaleoclimate.txt", sep="\t", row.names=FALSE, col.names=TRUE)


#############################################################################


# make regional and continental datasets from the gridded ones
# and add the information into Data_Regions and Data_Continents
# NOTE that this code overwrites the Data_Regions and Data_Continents files!

setwd("D://Data//Gammadiv//FilesForUpload")
se <- function(x) { # function for standard error of mean
  x <- x[!is.na(x)]
  return(sqrt(var(x)/length(x)))
}
lengthnotna <- function(x) (length(x[which(!is.na(x))])) # function to count the number of non-NA entries in a vector
luni <- function(x) (length(unique(x))) # function to count unique entries

clim <- read.delim("Data_GriddedPaleoclimate.txt", as.is=TRUE)
head(clim)
# sanity check - each cell only in one region?
regcon <- clim[,c(1:3)]
regcon <- regcon[which(!duplicated(regcon)),]
any(duplicated(regcon$WM_ID))
# get the variance column
clim$NPP_50 <- clim$NPP_50high - clim$NPP_50low
summary(clim$NPP_50)
hist(1/clim$NPP_50)
head(clim)
summary(as.factor(clim$AutomaticStage))
summary(as.factor(clim$BestStage))

# regional dataset

unireg <- unique(clim$Region) # the regions
unigltime <- unique(clim$BestStage[!is.na(clim$BestStage)]) # the time bins / stratigraphic stages
mydata <- c("AutoPlants", "BestPlants") # the types of dataset
clim$BestStage <- as.factor(clim$BestStage)
clim$AutomaticStage <- as.factor(clim$AutomaticStage)
for (i in seq(along=unireg)) { # loop through the regions
  dsi <- clim[which(clim$Region==unireg[i]),] # dataset for region i
  for (k in seq(along=mydata)) { # loop through the data types
    if (k==1) { # automatic stages
      dsik <- dsi[which(!is.na(dsi$AutomaticStage)),]
      names(dsik)[4] <- "MyStage"
    } else { # best stages
      dsik <- dsi[which(!is.na(dsi$BestStage)),]
      names(dsik)[5] <- "MyStage"
    }
    # only continue if the dataset contains rows
    if (length(dsik[,1])>0) {
      
      # get NPP summaries for each time bin 
      numloc <- tapply(dsik$NPP_hdr, dsik$MyStage, lengthnotna) # total number of grid cells with NPP data (ie both MAT and MAP) within each global time bin in region i
      meanhdr <- tapply(dsik$NPP_hdr, dsik$MyStage, mean, na.rm=TRUE) # mean NPP (from MAT and MAP hdr values within grid cells) within each global time bin in region i
      sehdr <- tapply(dsik$NPP_hdr, dsik$MyStage, se) # standard error of the mean NPP within each global time bin and in region i
      mean50low <- tapply(dsik$NPP_50low, dsik$MyStage, mean, na.rm=TRUE) # mean lower bound of the 50% interval for NPP (mean of the lower bounds within grid cells) within each global time bin in region i
      mean50high <- tapply(dsik$NPP_50high, dsik$MyStage, mean, na.rm=TRUE)
      meanmin <- tapply(dsik$NPP_meanmin, dsik$MyStage, mean, na.rm=TRUE) # mean minimum NPP (from the mean of the MAT and MAP minimum estimates within grid cells) within each global time bin in region i
      meanmax <- tapply(dsik$NPP_meanmax, dsik$MyStage, mean, na.rm=TRUE) # mean maximum NPP (from the mean of the MAT and MAP maximum estimates within grid cells) within each global time bin in region i
      unibin <- unique(as.character(dsik$MyStage))
      # have to loop through the bins for weighted means (bugger)
      for (j in seq(along=unibin)) {
        dsij <- dsik[which(dsik$MyStage==unibin[j]),] # dataset for region i, time bin k
        meanhdrwght50k <- weighted.mean(dsij$NPP_hdr, w=1/dsij$NPP_50, na.rm=TRUE) 
        meanminwghtk <- weighted.mean(dsij$NPP_meanmin, w=1/dsij$NPP_50, na.rm=TRUE) # mean of the minimum estimates, weighted by the inverse of the 50% density interval to get values comparable to hdr
        meanmaxwghtk <- weighted.mean(dsij$NPP_meanmax, w=1/dsij$NPP_50, na.rm=TRUE) # mean of the maximum estimates, weighted by the inverse of the 50% density interval to get values comparable to hdr
        names(meanhdrwght50k) <- unibin[j]
        names(meanminwghtk) <- unibin[j]
        names(meanmaxwghtk) <- unibin[j]
        if (j==1) {
          meanhdrwght50 <- meanhdrwght50k
          meanminwght50 <- meanminwghtk
          meanmaxwght50 <- meanmaxwghtk
        } else {
          meanhdrwght50 <-append(meanhdrwght50, meanhdrwght50k)
          meanminwght50 <- append(meanminwght50, meanminwghtk)
          meanmaxwght50 <- append(meanmaxwght50, meanmaxwghtk)
        }
      } # close loop through time bins
      # now pull all data together for NPP
      newlinen <- data.frame(Region=rep(unireg[i], length(unigltime)), StratigraphicStage=unigltime, NumCells=numloc[match(unigltime, names(numloc))], meanhdr=meanhdr[match(unigltime, names(meanhdr))], sehdr=sehdr[match(unigltime, names(sehdr))], meanhdrwght50=meanhdrwght50[match(unigltime, names(meanhdrwght50))], mean50low=mean50low[match(unigltime, names(mean50low))], mean50high=mean50high[match(unigltime, names(mean50high))], meanmin=meanmin[match(unigltime, names(meanmin))], meanminwght50=meanminwght50[match(unigltime, names(meanminwght50))], meanmax=meanmax[match(unigltime, names(meanmax))], meanmaxwght50=meanmaxwght50[match(unigltime, names(meanmaxwght50))])
      names(newlinen)[3:length(newlinen)] <- paste("NPP", mydata[k], names(newlinen)[3:length(newlinen)], sep="_")
    } else { # if the dataset contains no rows - no NPP data at all for a region!
      newlinen <- data.frame(Region=rep(unireg[i], length(unigltime)), StratigraphicStage=unigltime, NumCells=0, meanhdr=NA, sehdr=NA, meanhdrwght50=NA, mean50low=NA, mean50high=NA, meanmin=NA, meanminwght50=NA, meanmax=NA, meanmaxwght50=NA)
      names(newlinen)[3:length(newlinen)] <- paste("NPP", mydata[k], names(newlinen)[3:length(newlinen)], sep="_")
    }
    # unite the newlinen datasets across all data types
    if (k==1) (newline <- newlinen) else (newline <- merge(newline, newlinen, all=TRUE))
  } # loop through the data types
  # then merge all the region data together
  if (i==1) (ds <- newline) else (ds <- rbind(ds, newline))
}

summary(ds)
unique(ds$Region)
ds[which(ds$NPP_BestPlants_NumCells==0),]
ds[which(is.na(ds$NPP_BestPlants_NumCells)),]
# sort out some NAs
names(ds)
cellcols <- grep("NumCells", names(ds))
for (i in seq(along=cellcols)) (ds[which(is.na(ds[,cellcols[i]])),cellcols[i]] <- 0)
summary(ds)
othercols <- c(3:length(ds))
othercols <- othercols[!is.element(othercols, cellcols)]
for (i in seq(along=othercols)) (ds[which(is.nan(ds[,othercols[i]])),othercols[i]] <- NA)
table(ds$Region, ds$StratigraphicStage)

# unite with the mammalian dataset
dsr <- read.delim("Data_Regions.txt", as.is=TRUE)
head(dsr)
# note that this code works on the assumption that Data_Regions is the dataset 
# generated with Rscript_MammalDataProcessing.r. If you work with the uploaded
# dataset, use the next line to reset the data.
# dsr <- dsr[,1:25] # gets rid of NPP columns if these have already been added
summary(as.factor(dsr$Region))
summary(as.factor(dsr$GlobalTimeBins))
dsr2 <- merge(dsr, ds, by=c(1:2), all=TRUE, sort=FALSE)
table(dsr2$Region, dsr2$GlobalTimeBins)
summary(as.factor(dsr2$Region))
summary(as.factor(dsr2$GlobalTimeBins))
head(dsr2)
# IF YOU DO NOT WANT TO OVERWRITE THE FILE UPLOADED WITH THE PAPER,
# change the file name below.
write.table(dsr2, "Data_Regions.txt", sep="\t", row.names=FALSE, col.names=TRUE)


# continental dataset

unireg <- unique(clim$Continent)
for (i in seq(along=unireg)) { # loop through the regions
  dsi <- clim[which(clim$Continent==unireg[i]),] # dataset for continent i
  for (k in seq(along=mydata)) { # loop through the data types
    if (k==1) { # automatic stages
      dsik <- dsi[which(!is.na(dsi$AutomaticStage)),]
      names(dsik)[4] <- "MyStage"
    } else { # best stages
      dsik <- dsi[which(!is.na(dsi$BestStage)),]
      names(dsik)[5] <- "MyStage"
    }
    # only continue if the dataset contains rows
    if (length(dsik[,1])>0) {
      
      # get NPP summaries for each time bin 
      numloc <- tapply(dsik$NPP_hdr, dsik$MyStage, lengthnotna) # total number of grid cells with NPP data (ie both MAT and MAP) within each global time bin in region i
      meanhdr <- tapply(dsik$NPP_hdr, dsik$MyStage, mean, na.rm=TRUE) # mean NPP (from MAT and MAP hdr values within grid cells) within each global time bin in region i
      sehdr <- tapply(dsik$NPP_hdr, dsik$MyStage, se) # standard error of the mean NPP within each global time bin and in region i
      mean50low <- tapply(dsik$NPP_50low, dsik$MyStage, mean, na.rm=TRUE) # mean lower bound of the 50% interval for NPP (mean of the lower bounds within grid cells) within each global time bin in region i
      mean50high <- tapply(dsik$NPP_50high, dsik$MyStage, mean, na.rm=TRUE)
      meanmin <- tapply(dsik$NPP_meanmin, dsik$MyStage, mean, na.rm=TRUE) # mean minimum NPP (from the mean of the MAT and MAP minimum estimates within grid cells) within each global time bin in region i
      meanmax <- tapply(dsik$NPP_meanmax, dsik$MyStage, mean, na.rm=TRUE) # mean maximum NPP (from the mean of the MAT and MAP maximum estimates within grid cells) within each global time bin in region i
      unibin <- unique(as.character(dsik$MyStage))
      # have to loop through the bins for weighted means (bugger)
      for (j in seq(along=unibin)) {
        dsij <- dsik[which(dsik$MyStage==unibin[j]),] # dataset for region i, time bin k
        meanhdrwght50k <- weighted.mean(dsij$NPP_hdr, w=1/dsij$NPP_50, na.rm=TRUE) 
        meanminwghtk <- weighted.mean(dsij$NPP_meanmin, w=1/dsij$NPP_50, na.rm=TRUE) # mean of the minimum estimates, weighted by the inverse of the 50% density interval to get values comparable to hdr
        meanmaxwghtk <- weighted.mean(dsij$NPP_meanmax, w=1/dsij$NPP_50, na.rm=TRUE) # mean of the maximum estimates, weighted by the inverse of the 50% density interval to get values comparable to hdr
        names(meanhdrwght50k) <- unibin[j]
        names(meanminwghtk) <- unibin[j]
        names(meanmaxwghtk) <- unibin[j]
        if (j==1) {
          meanhdrwght50 <- meanhdrwght50k
          meanminwght50 <- meanminwghtk
          meanmaxwght50 <- meanmaxwghtk
        } else {
          meanhdrwght50 <-append(meanhdrwght50, meanhdrwght50k)
          meanminwght50 <- append(meanminwght50, meanminwghtk)
          meanmaxwght50 <- append(meanmaxwght50, meanmaxwghtk)
        }
      } # close loop through time bins
      # now pull all data together for NPP
      newlinen <- data.frame(Region=rep(unireg[i], length(unigltime)), StratigraphicStage=unigltime, NumCells=numloc[match(unigltime, names(numloc))], meanhdr=meanhdr[match(unigltime, names(meanhdr))], sehdr=sehdr[match(unigltime, names(sehdr))], meanhdrwght50=meanhdrwght50[match(unigltime, names(meanhdrwght50))], mean50low=mean50low[match(unigltime, names(mean50low))], mean50high=mean50high[match(unigltime, names(mean50high))], meanmin=meanmin[match(unigltime, names(meanmin))], meanminwght50=meanminwght50[match(unigltime, names(meanminwght50))], meanmax=meanmax[match(unigltime, names(meanmax))], meanmaxwght50=meanmaxwght50[match(unigltime, names(meanmaxwght50))])
      names(newlinen)[3:length(newlinen)] <- paste("NPP", mydata[k], names(newlinen)[3:length(newlinen)], sep="_")
    } else { # if the dataset contains no rows - no NPP data at all for a region!
      newlinen <- data.frame(Region=rep(unireg[i], length(unigltime)), StratigraphicStage=unigltime, NumCells=0, meanhdr=NA, sehdr=NA, meanhdrwght50=NA, mean50low=NA, mean50high=NA, meanmin=NA, meanminwght50=NA, meanmax=NA, meanmaxwght50=NA)
      names(newlinen)[3:length(newlinen)] <- paste("NPP", mydata[k], names(newlinen)[3:length(newlinen)], sep="_")
    }
    # unite the newlinen datasets across all data types
    if(k==1) (newline <- newlinen) else (newline <- merge(newline, newlinen, all=TRUE))
  } # loop through the data types
  # then merge all the region data together
  if (i==1) (ds <- newline) else (ds <- rbind(ds, newline))
}

summary(ds)
names(ds)[1] <- "Continent"
unique(ds$Continent)
ds[which(ds$NPP_AutoPlants_NumCells==0),]
ds[which(is.na(ds$NPP_AutoPlants_NumCells)),]
# sort out some NAs
names(ds)
cellcols <- grep("NumCells", names(ds))
for (i in seq(along=cellcols)) (ds[which(is.na(ds[,cellcols[i]])),cellcols[i]] <- 0)
summary(ds)
othercols <- c(3:length(ds))
othercols <- othercols[!is.element(othercols, cellcols)]
for (i in seq(along=othercols)) (ds[which(is.nan(ds[,othercols[i]])),othercols[i]] <- NA)
table(ds$Continent, ds$StratigraphicStage)

# unite with the mammalian dataset
dsc <- read.delim("Data_Continents.txt", as.is=TRUE)
head(dsc)
# note that this code works on the assumption that Data_Continents is the dataset 
# generated with Rscript_MammalDataProcessing.r. If you work with the uploaded
# dataset, use the next line to reset the data.
# dsc <- dsc[,1:24] # gets rid of NPP columns if these have already been added
summary(as.factor(dsc$Continent))
summary(as.factor(dsc$GlobalTimeBins))
dsc2 <- merge(dsc, ds, by=c(1:2), all=TRUE, sort=FALSE)
table(dsc2$Continent, dsc2$GlobalTimeBins)
summary(as.factor(dsc2$Continent))
summary(as.factor(dsc2$GlobalTimeBins))
head(dsc2)
# IF YOU DO NOT WANT TO OVERWRITE THE FILE UPLOADED WITH THE PAPER,
# change the file name below.
write.table(dsc2, "Data_Continents.txt", sep="\t", row.names=FALSE, col.names=TRUE)

