Monthly climate data at 1 km resolution come from down-scaling terraClimate to 1 km with CHELSA and projecting CHELSA forward and backward using terraClimate. This file is here: clark.unix/getClimateData.

This code runs from clark.unix. Replace my path with yours to the clark.unix directory:

user <- 'mac2unix' 
if(user == 'mac2unix'){
  path <- "/Volumes/research/clark/clark.unix/"
}else if(user == 'mac'){
  path <- "/Users/jimclark/makeMastOnJimClark/"
}else if(user == 'mac64'){
  path <- "~/Library/Mobile Documents/com~apple~CloudDocs/Documents/makeMastOnJimClark/"
}

Set paths, source functions:

apath <- paste( path, 'makeMast/', sep='' )
hpath <- paste( apath, 'climateBuild/', sep='' )
tpath <- paste( hpath, 'terraClimateFiles/', sep='' )

source( paste(hpath, 'climateFunctions.r', sep='') )
source( paste(apath, 'RFunctions/mastifFunctions.r', sep='') )

Enter the locations (lon, lat) in the matrix xy, give them rownames, and specify vars (some combination of tmin, tmax, ppt):

xy <- matrix( c( -79.0944674, 35.98487734, -110.45039, 44.96276, 31.5969, -24.9948), 3, 2, byrow=T)
rownames(xy) <- c( 'Blackwood', 'Lamar Valley', 'Skukuza' )
colnames(xy) <- c('lon','lat')
vars <- c("tmin","tmax","ppt")

Here they are on a map:

Map of locations in Yellowstone (Lamar Valley), Duke Forest (Blackwood), and Kruger (Skukuza).

Map of locations in Yellowstone (Lamar Valley), Duke Forest (Blackwood), and Kruger (Skukuza).

Here is the xy matrix:

Specify the years, months, and tpath in the call to getTileTerraClimate. Here I get all months for the last 10 years. The most recent year is now 2019:

library( terra )
startYr <- 1990
endYr   <- 2021
tmp <- getTileTerraClimate(lonLat = xy, years = startYr:endYr, months = 1:12, vars = vars, 
                           DOWNSCALE = T, path = tpath) # always set DOWNSCALE = T
tall <- tmp$tallFormat
wide <- tmp$wideFormat

There are two formats. The tall format has all variables in vars. The wide format is a list, one for each element of vars.

Note that the wide format has location in the first two columns. Units for ppt are mm. You can get approximate temp like this:

temp <- (wide$tmin + wide$tmax)/2

If you want PET, first get ppt and temp. Then do this:

prec <- wide$ppt

# remove lon, lat from first two columns
temp <- temp[,-c(1:2)]
prec <- prec[,-c(1:2)]

# year and month indices for columns
cc <- columnSplit(colnames(temp),'_')
yi <- as.numeric(cc[,1])
mi <- as.numeric(cc[,2])

lat <- xy[,'lat']  # latitude needed for daylength
pet <- monthlyPET(yi, mi, temp, prec, lat=lat)[[1]]   
rownames(pet) <- rownames(temp)
colnames(pet) <- colnames(temp)
def  <- pet - prec

Here are some quick plots, check data for your area.

Moisture deficits (orange) integrate changes in temperature (PET) and precipitation.

Moisture deficits (orange) integrate changes in temperature (PET) and precipitation.

I’ve used blue to remind that a negative deficit is ‘wet’.

Below is code that can be used to draw the water balance as norms for a given location. Here I evaluate monthly water balance for two eight-yr intervals followed by the difference between those intervals.

par( mfrow = c(3, 3), bty = 'n', mar = c(2,2,1,1), omi = c(.5, .5, .1, .5) )

c1 <- which( yr <= 1997 )
c2 <- which( yr >= 2014 )
cc <- c1
ik <- 0

for( i in 1:6){
  
  ik <- ik + 1
  
  if( i == 4)ik <- 1
  if( i >= 4 )cc <- c2
  pe <- tapply( pet[ik,cc], mo[cc], mean )
  ti <- tapply( temp[ik,cc], mo[cc], mean )
  pr <- tapply( prec[ik,cc], mo[cc], mean )
  de <- tapply( def[ik,mo[cc]], mo[cc], mean )
  
  plot( 1:12, pr, type = 'l', ylim = c(-20, 170), xlab = '', ylab = '', col = 'darkblue', lwd = 2 )
  lines( 1:12, ti, col = 'red', lwd = 2 )
  lines( 1:12, pe, col = '#800026', lwd = 2 )
  
  # shade polygons
  pts <- intersectLines( 1:12, pe, pr )
  if( length(pts$poly1) > 0 ){
    for( j in 1:length( pts$poly1 ) )
      polygon(pts$poly1[[j]][,1], pts$poly1[[j]][,2], col = '#fed976')
  }
  if( length(pts$poly2) > 0 ){
    for( j in 1:length( pts$poly2 ) )
      polygon(pts$poly2[[j]][,1], pts$poly2[[j]][,2], col = '#9ecae1')
  }
  segments( 1:12, pe, 1:12, pr, col = '#bd0026' )
  
  abline( h = 0, lty = 2, col = 'grey', lwd = 2 )
  title( rownames(pet)[i] )
  
  if( i == 3 )mtext( '1990 - 1997', 4, line = 1 )
  if( i == 1 ){
    text( 5, ti[5] + 15, 'T', col = 'red', cex = 1.4 )
    text( 2, pr[2] + 25, 'P', col = 'darkblue', cex = 1.4 )
    text( 5, pe[5] + 35, 'PET', col = '#800026', cex = 1.4 )
  }
  if( i == 6 )mtext( '2014 - 2021', 4, line = 1 )
}

for( i in 1:3 ){
  de <- tapply( def[i,c2], mo[c2], mean ) - tapply( def[i,c1], mo[c1], mean )
  plot( 1:12, de, ylim = c( -30, 30), type = 'l', lwd = 2 )
  abline( h = 0, lty = 2, col = 'grey', lwd = 2 )
  
  pts <- intersectLines( 1:12, de, de*0 )
  if( length(pts$poly1) > 0 ){
  for( j in 1:length( pts$poly1 ) )
    polygon(pts$poly1[[j]][,1], pts$poly1[[j]][,2], col = '#fed976')
  }
  if( length(pts$poly2) > 0 ){
  for( j in 1:length( pts$poly2 ) )
    polygon(pts$poly2[[j]][,1], pts$poly2[[j]][,2], col = '#9ecae1')
  }
  segments( 1:12, de*0, 1:12, de, col = '#bd0026' )
}
mtext( 'Difference', 4, line = 1, outer = F)
mtext( 'Month', 1, line = 1, outer = T)
mtext( 'mm', 2, line = 1, outer = T )
Water balance shows deficits where P < PET (orange) and surpluses when P > PET (blue). Rows show norms for two 8-yr periods. The difference plots at bottom show monthly change, with increased deficit as positive values.

Water balance shows deficits where P < PET (orange) and surpluses when P > PET (blue). Rows show norms for two 8-yr periods. The difference plots at bottom show monthly change, with increased deficit as positive values.