R Code for Climate Mathematics:

Theory and Applications

 

A Cambridge University Press book By

SSP Shen and RCJ Somerville

 

 

Version 1.0 released in July 2019 and coded by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
https://shen.sdsu.edu
Email:

 

Version 2.0 compiled and updated by Momtaza Sayd
San Diego State University May 2021
Version 3.0 compiled and updated by Joaquin Stawsky
San Diego State University June 2022

 

 

 

Chapter 2: Basics of R Programming

R as a smart calculator

1+4
## [1] 5
2+pi/4-0.8
## [1] 1.985398
x <- 1
y <- 2
z <- 4
t <- 2*x^y-z
t
## [1] -2
u = 2        # "=" sign and "<-" are almost equivalent
v = 3        # The text behind the "#" sign is comments
u+v
## [1] 5
sin(u*v)    # u*v = 6 in the sine function is considered a radian by R 
## [1] -0.2794155

 

Define a sequence in R

#Enter temperature data in c() 
tmax <- c(77, 72, 75, 73, 66, 64, 59)
#Show the data
tmax
## [1] 77 72 75 73 66 64 59
#Generate same sequence using different methods
seq(1,8)
## [1] 1 2 3 4 5 6 7 8
seq(8)
## [1] 1 2 3 4 5 6 7 8
seq(1,8, by = 1)
## [1] 1 2 3 4 5 6 7 8
seq(1,8, length = 8)
## [1] 1 2 3 4 5 6 7 8
seq(1,8, length.out = 8)
## [1] 1 2 3 4 5 6 7 8

 

Define a function in R

#Define a function
samfctn <- function(x) x*x
samfctn(4)
## [1] 16
fctn2 <- function(x,y,z) x+y-z/2
fctn2(1,2,3)
## [1] 1.5

 

Plot with R

#Plot temperature data
plot(1:7, c(77, 72, 75, 73, 66, 64, 59))

#More plot examples
plot(sin, -pi, 2*pi)   #plot the curve of y = sin(x) from -pi to 2 pi

square <- function(x) x*x   #Define a function
plot(square, -3, 2)   # Plot the defined function

# Plot a 3D surface
x <- seq(-1, 1, length = 100)
y <- seq(-1, 1, length = 100)
z <- outer(x, y, function(x, y)(1-x^2-y^2))  
# 'outer(x,y, function)' renders z function on the x, y grid
persp(x,y,z, theta = 330) 

# yields a 3D surface with perspective angle 330 deg

#Contour plot
contour(x,y,z) #lined contours

filled.contour(x,y,z) #color map of contours with default colors

 

Symbolic calculations by R

D(expression(x^2,'x'), 'x') 
## 2 * x
# Take derivative of x^2 w.r.t. x 
2 * x #The answer is 2x
##   [1] -2.00000000 -1.95959596 -1.91919192 -1.87878788 -1.83838384 -1.79797980
##   [7] -1.75757576 -1.71717172 -1.67676768 -1.63636364 -1.59595960 -1.55555556
##  [13] -1.51515152 -1.47474747 -1.43434343 -1.39393939 -1.35353535 -1.31313131
##  [19] -1.27272727 -1.23232323 -1.19191919 -1.15151515 -1.11111111 -1.07070707
##  [25] -1.03030303 -0.98989899 -0.94949495 -0.90909091 -0.86868687 -0.82828283
##  [31] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263 -0.58585859
##  [37] -0.54545455 -0.50505051 -0.46464646 -0.42424242 -0.38383838 -0.34343434
##  [43] -0.30303030 -0.26262626 -0.22222222 -0.18181818 -0.14141414 -0.10101010
##  [49] -0.06060606 -0.02020202  0.02020202  0.06060606  0.10101010  0.14141414
##  [55]  0.18181818  0.22222222  0.26262626  0.30303030  0.34343434  0.38383838
##  [61]  0.42424242  0.46464646  0.50505051  0.54545455  0.58585859  0.62626263
##  [67]  0.66666667  0.70707071  0.74747475  0.78787879  0.82828283  0.86868687
##  [73]  0.90909091  0.94949495  0.98989899  1.03030303  1.07070707  1.11111111
##  [79]  1.15151515  1.19191919  1.23232323  1.27272727  1.31313131  1.35353535
##  [85]  1.39393939  1.43434343  1.47474747  1.51515152  1.55555556  1.59595960
##  [91]  1.63636364  1.67676768  1.71717172  1.75757576  1.79797980  1.83838384
##  [97]  1.87878788  1.91919192  1.95959596  2.00000000
fx <- expression(x^2,'x')  #assign a function  
D(fx,'x') #differentiate the function w.r.t. x
## 2 * x
2 * x  #The answer is 2x
##   [1] -2.00000000 -1.95959596 -1.91919192 -1.87878788 -1.83838384 -1.79797980
##   [7] -1.75757576 -1.71717172 -1.67676768 -1.63636364 -1.59595960 -1.55555556
##  [13] -1.51515152 -1.47474747 -1.43434343 -1.39393939 -1.35353535 -1.31313131
##  [19] -1.27272727 -1.23232323 -1.19191919 -1.15151515 -1.11111111 -1.07070707
##  [25] -1.03030303 -0.98989899 -0.94949495 -0.90909091 -0.86868687 -0.82828283
##  [31] -0.78787879 -0.74747475 -0.70707071 -0.66666667 -0.62626263 -0.58585859
##  [37] -0.54545455 -0.50505051 -0.46464646 -0.42424242 -0.38383838 -0.34343434
##  [43] -0.30303030 -0.26262626 -0.22222222 -0.18181818 -0.14141414 -0.10101010
##  [49] -0.06060606 -0.02020202  0.02020202  0.06060606  0.10101010  0.14141414
##  [55]  0.18181818  0.22222222  0.26262626  0.30303030  0.34343434  0.38383838
##  [61]  0.42424242  0.46464646  0.50505051  0.54545455  0.58585859  0.62626263
##  [67]  0.66666667  0.70707071  0.74747475  0.78787879  0.82828283  0.86868687
##  [73]  0.90909091  0.94949495  0.98989899  1.03030303  1.07070707  1.11111111
##  [79]  1.15151515  1.19191919  1.23232323  1.27272727  1.31313131  1.35353535
##  [85]  1.39393939  1.43434343  1.47474747  1.51515152  1.55555556  1.59595960
##  [91]  1.63636364  1.67676768  1.71717172  1.75757576  1.79797980  1.83838384
##  [97]  1.87878788  1.91919192  1.95959596  2.00000000
fx <- expression(x^2*sin(x),'x') 
#Change the expression and use the same derivative command
D(fx,'x')
## 2 * x * sin(x) + x^2 * cos(x)
2 * x * sin(x) + x^2 * cos(x)
##   [1] 2.2232442755 2.1621237197 2.1001569223 2.0374552454 1.9741295366
##   [6] 1.9102900272 1.8460462307 1.7815068435 1.7167796456 1.6519714028
##  [11] 1.5871877703 1.5225331975 1.4581108337 1.3940224359 1.3303682779
##  [16] 1.2672470605 1.2047558237 1.1429898609 1.0820426339 1.0220056908
##  [21] 0.9629685847 0.9050187953 0.8482416517 0.7927202577 0.7385354186
##  [26] 0.6857655712 0.6344867148 0.5847723450 0.5366933900 0.4903181485
##  [31] 0.4457122306 0.4029385012 0.3620570247 0.3231250141 0.2861967807
##  [36] 0.2513236874 0.2185541047 0.1879333686 0.1595037417 0.1333043769
##  [41] 0.1093712835 0.0877372965 0.0684320482 0.0514819428 0.0369101336
##  [46] 0.0247365036 0.0149776477 0.0076468594 0.0027541183 0.0003060825
##  [51] 0.0003060825 0.0027541183 0.0076468594 0.0149776477 0.0247365036
##  [56] 0.0369101336 0.0514819428 0.0684320482 0.0877372965 0.1093712835
##  [61] 0.1333043769 0.1595037417 0.1879333686 0.2185541047 0.2513236874
##  [66] 0.2861967807 0.3231250141 0.3620570247 0.4029385012 0.4457122306
##  [71] 0.4903181485 0.5366933900 0.5847723450 0.6344867148 0.6857655712
##  [76] 0.7385354186 0.7927202577 0.8482416517 0.9050187953 0.9629685847
##  [81] 1.0220056908 1.0820426339 1.1429898609 1.2047558237 1.2672470605
##  [86] 1.3303682779 1.3940224359 1.4581108337 1.5225331975 1.5871877703
##  [91] 1.6519714028 1.7167796456 1.7815068435 1.8460462307 1.9102900272
##  [96] 1.9741295366 2.0374552454 2.1001569223 2.1621237197 2.2232442755
fxy <- expression(x^2+y^2, 'x','y') 
#One can define a function of 2 or more variables
fxy #Renders an expression of the function in terms of x and y
## expression(x^2 + y^2, "x", "y")
#expression(x^2 + y^2, "x", "y")
D(fxy,'x') #yields the partial derivative with respect to x: 2 * x
## 2 * x
D(fxy,'y') #yields the partial derivative with respect to y: 2 * y
## 2 * y
square <- function(x) x^2
#Integrate x^2 from 0 to 1 equals to 1/3 with details below
integrate(square,0,1)
## 0.3333333 with absolute error < 3.7e-15
#0.3333333 with absolute error < 3.7e-15

#Integrate cos(x) from 0 to pi/2 equals to 1 with details below
integrate(cos,0,pi/2) 
## 1 with absolute error < 1.1e-14
#1 with absolute error < 1.1e-14

 

 

Simple statistics by R

x <- rnorm(10) #generate 10 normally distributed numbers
x
##  [1]  1.23263330  0.01294500 -0.95786978 -0.47574866 -0.33991132 -0.39052728
##  [7]  0.78613457  0.08218826 -2.55529275 -0.12642525
mean(x)
## [1] -0.2731874
var(x)
## [1] 1.03949
sd(x)
## [1] 1.019554
median(x)
## [1] -0.2331683
quantile(x)
##          0%         25%         50%         75%        100% 
## -2.55529275 -0.45444332 -0.23316829  0.06487744  1.23263330
range(x) #yields the min and max of x
## [1] -2.555293  1.232633
max(x)
## [1] 1.232633
boxplot(x) #yields the box plot of x

w <- rnorm(1000)
hist(w) 

#yields the histogram of 1000 random numbers with a normal distribution

summary(rnorm(12)) #statistical summary of the data sequence
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -1.0572 -0.2931  0.2651  0.1315  0.6012  1.0083
#Linear regression and linear trend line
#2007-2016 data of the global temperature anomalies
#Source: NOAAGlobalTemp data
t <- 2007:2016
T <- c(.36,.30, .39, .46, .33, .38, .42, .50, .66, .70)
lm(T ~ t) #Linear regression model of temp vs time
## 
## Call:
## lm(formula = T ~ t)
## 
## Coefficients:
## (Intercept)            t  
##   -73.42691      0.03673
#Temperature change rate is 0.03673  deg C/yr or 0.37 deg C/decade
# expression(paste(...)) is one way of using the degree symbol
# ...pasting the ° symbol directly is another way
plot(t,T, type = "o",xlab = "Year",ylab = expression(paste("Temperature [",~degree,"C]")), 
     main = expression(paste("2007-2016 Global Temperature Anomalies\n and Their Linear Trend [0.37 °C/decade]")))
abline(lm(T ~ t), lwd = 2, col = "red") #Regression line

 

YouTube tutorial: Input data by reading a csv file into R

#The R packages and the datasets used in this book are 
#listed below and can be downloaded and installed first 
#before proceeding to the R codes in the rest of the book. 
#The R packages: 
#animation, chron, e1071, fields, ggplot2, lattice, 
#latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf, 
#NLRoot, RColorBrewer, rgdal, rasterVis, raster, sp, TTR

#To load the package "animation", you can do 
library(animation)

#You can also load all these packages in one shot
#using pacman
#install.packages("pacman")
library(pacman)
pacman::p_load(animation, chron, e1071, fields, ggplot2, lattice, 
            latticeExtra, maps, mapdata, mapproj, matrixStats, ncdf4, 
            NLRoot, RColorBrewer, rgdal, rasterVis, raster, sp, TTR)
## also installing the dependencies 'dotCall64', 'spam'
## 
## The downloaded binary packages are in
##  /var/folders/q8/xzb0sp_n4w12lbf3bfc2vf900000gp/T//RtmpEfo0WW/downloaded_packages
## 
## fields installed
## 
## The downloaded binary packages are in
##  /var/folders/q8/xzb0sp_n4w12lbf3bfc2vf900000gp/T//RtmpEfo0WW/downloaded_packages
## 
## mapdata installed
## 
## The downloaded binary packages are in
##  /var/folders/q8/xzb0sp_n4w12lbf3bfc2vf900000gp/T//RtmpEfo0WW/downloaded_packages
## 
## mapproj installed
## 
## The downloaded binary packages are in
##  /var/folders/q8/xzb0sp_n4w12lbf3bfc2vf900000gp/T//RtmpEfo0WW/downloaded_packages
## 
## NLRoot installed
## 
## The downloaded binary packages are in
##  /var/folders/q8/xzb0sp_n4w12lbf3bfc2vf900000gp/T//RtmpEfo0WW/downloaded_packages
## 
## rgdal installed
## also installing the dependency 'hexbin'
## 
## The downloaded binary packages are in
##  /var/folders/q8/xzb0sp_n4w12lbf3bfc2vf900000gp/T//RtmpEfo0WW/downloaded_packages
## 
## rasterVis installed
#The zipped data file:
# https://www.cambridge.org/climatemathematics/data.zip

#On your computer, you can create a directory called 
#climmath under your user name. 
#The one used in the book is Users/sshen/climmath  
#You unzip the data and move the data folder under 
#the Users/sshen/climmath directory. 
#A data folder will created:  
#Users/sshen/climmath/data. 
#The data folder contains about 400 MB of data. 
#Place all the R codes in the directory  Users/sshen/climmath. 
#Then, you can run all the codes in this book after replacing sshen 
#by your user name on your own computer.