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
Email: sshen@sdsu.edu
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] 52+pi/4-0.8## [1] 1.985398x <- 1
y <- 2
z <- 4
t <- 2*x^y-z
t## [1] -2u = 2 # "=" sign and "<-" are almost equivalent
v = 3 # The text behind the "#" sign is comments
u+v## [1] 5sin(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 8seq(8)## [1] 1 2 3 4 5 6 7 8seq(1,8, by = 1)## [1] 1 2 3 4 5 6 7 8seq(1,8, length = 8)## [1] 1 2 3 4 5 6 7 8seq(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] 16fctn2 <- 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.00000000fx <- expression(x^2,'x') #assign a function
D(fx,'x') #differentiate the function w.r.t. x## 2 * x2 * 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.00000000fx <- 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.2232442755fxy <- 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 * xD(fxy,'y') #yields the partial derivative with respect to y: 2 * y## 2 * ysquare <- 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] -0.40518274 0.43551556 -1.19994395 -0.85516984 -0.02758804 -0.04578632
## [7] -0.16097080 -0.33794435 -1.16459046 -0.34227721mean(x)## [1] -0.4103938var(x)## [1] 0.2730219sd(x)## [1] 0.522515median(x)## [1] -0.3401108quantile(x)## 0% 25% 50% 75% 100%
## -1.19994395 -0.74267306 -0.34011078 -0.07458244 0.43551556range(x) #yields the min and max of x## [1] -1.1999440 0.4355156max(x)## [1] 0.4355156boxplot(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.93697 -0.44464 -0.18273 -0.01719 -0.02869 2.17169#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)
#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.