1.1 Why use R?

R was initially developed by Robert Gentleman and Ross Ihaka of the Department of Statistics at the University of Auckland. R is increasingly becoming the default software package in many areas of science. There are a number of reasons for this:

• R includes a very large number of tools, functions and packages
  • packages are libraries of tools
  • packages are written by scientists in different subject areas
• R has the latest methods and tools
• New tools are in R 10-20 years before commercial software
  • e.g. GWR was around for 10 years before being in ArcGIS
• The tools in R are open (i.e. the source code is visible)
  • you can see exactly what is being done
  • most commercial tools are hidden - black boxes
• R is free!

For these reasons R is becoming widely used in many areas of scientific activity and quantitative research.

R can be found at the CRAN website:

1.2 Working in R

There are 2 key points about working in R

1. When working in R, either writing your own code or copy and pasting from these materials, you should* write the code into a script or document. Go to File > New File > R Script** to open a new R file.

The reasons for this are so that you get used to using the R console and running the code will help your understanding of the code’s functionality. Then in order to run the code in the R console,a quick way to enter it is to highlight the code (with the mouse or using the keyboard controls) and the press ctrl-R or cmd-Enter on a Mac.

2. Learning is R is learning to drive. You may pass your test but ti become a good driver it is time behind the wheel that counts. The importance of learning by doing and getting your hands dirty cannot be overstated. Some of the code might look a bit fearsome when first viewed, especially in later session BUT the only really effective way to understand it is to give it a try. A further minor point is that in the code comments are prefixed by \(#\) and are ignored by R when entered into the console.

There are some core concepts that you should be familiar with:

1. Assignment: this is the basic process of giving R objects values:

vals <- c(4.3,7.1,6.3,5.2,3.2,2.1)

2. Operations: having assigned values to object that can be manipulated using functions:

vals*2

## [1]  8.6 14.2 12.6 10.4  6.4  4.2

sum(vals)

## [1] 28.2

mean(vals)

## [1] 4.7

3. Indexing: elements of R objects can be can be referred to:

vals[1]    # first element

## [1] 4.3

vals[1:3]   # a subset of elements 1 to 3

## [1] 4.3 7.1 6.3

sqrt(vals[1:3]) #square roots of the subset

## [1] 2.073644 2.664583 2.509980

vals[c(5,3,2)]  # a subset of elements 5,3,2 - note the ordering

## [1] 3.2 6.3 7.1

Note you can copy and paste the code above into a new R/RStudio script.

1.3 R packages

When you install R / RStudio it comes with a large number of tools already (refereed to as base functionality).

However, one of the joys of R, is the community of users. Users share what they do and create in R in a number of ways. One of these is through packages. Packages are collections of related functions that have been created, tested and supported with help files. These are bundled into a package and shared with other R users via the that users can download from the CRAN repository.

There are 1000s of packages in R. These contain set of tools and can be written by anyone. The number of packages is continually growing. When packages are installed these can be called as libraries. The background to R, along with documentation and information about packages as well as the contributors, can be found at the R Project website http://www.r-project.org.

Packages can be found at the CRAN website - https://cran.r-project.org/web/packages/:

Users install the package once to mount it on their computer, and then it can be called in R scripts as required.

The basic operations are

1. install the package before the first time it is used
   • you may have to set a mirror site the fist time you install a a package
   • this is only done once
2. load the package using the library function to use the package tools
   • this is done for each R session

So a typical way to do this is

install.packages("<package name>"", dep = TRUE)

For this workshop you will require a number of packages to be installed.

The code below will make sure you have all of the R packages you will need installed. The tests check whether a given package is in the list of installed packages, and installs it if not. Package installation only needs to be done once - not every time you use R - so if you are sure you have all of the packages installed already, you can skip this step.

# set a mirror
options(repos = c(CRAN = "http://cran.rstudio.com"))
# test for package existance and install
if (!is.element("tidyverse", installed.packages()))
    install.packages("tidyverse", dep = T)
if (!is.element("GISTools", installed.packages()))
    install.packages("GISTools", dep = T)
if (!is.element("GWmodel", installed.packages())) 
    install.packages("GWmodel", dep = T)
if (!is.element("MASS", installed.packages())) 
    install.packages("MASS", dep = T)
if (!is.element("rgdal", installed.packages()))
    install.packages("rgdal", dep = T)
if (!is.element("tmap", installed.packages())) 
    install.packages("tmap", dep = T)
if (!is.element("e1071", installed.packages()))
    install.packages("e1071", dep = T)
if (!is.element("knitr", installed.packages()))
    install.packages("knitr", dep = T)
if (!is.element("car", installed.packages()))
    install.packages("car", dep = T)
if (!is.element("gclus", installed.packages()))
    install.packages("gclus", dep = T)
if (!is.element("kableExtra", installed.packages()))
    install.packages("kableExtra", dep = T)

Note that the code above installs each package and all of is dependencies (dep = TRUE). Once we are sure all the packages are installed, you need to load them into the current session:

# load into the R session
library(tidyverse)
library(GISTools)
library(GWmodel)
library(MASS)
library(rgdal)
library(tmap)
library(e1071)
library(knitr)
library(car)
library(gclus)
library(kableExtra)

Also note that packages can be also be downloaded using the RStudio menu

Tools > Install Packages…

2.1 Loading tabular data from your your computer

## Mac
setwd("/Users/geoaco/Desktop/")
## Windows
setwd("C:\\")

dat <- read.csv("GWPCA_SA_2018_Workshopdata.csv")

## dimensions - rows and columns
dim(dat)

## [1] 1070   33

## column / variable names
names(dat)

##  [1] "ID"                   "Easting"              "Northing"            
##  [4] "Field_No"             "Field_East_Centroid"  "Field_North_Centroid"
##  [7] "Farmlet_No"           "Soil_Class_No"        "Elevation"           
## [10] "Slope"                "Aspect"               "Al"                  
## [13] "As"                   "Ca"                   "Cd"                  
## [16] "Co"                   "Cr"                   "Cu"                  
## [19] "Fe"                   "K"                    "Mg"                  
## [22] "Mn"                   "Mo"                   "Na"                  
## [25] "Ni"                   "P"                    "Pb"                  
## [28] "S"                    "Se"                   "Ti"                  
## [31] "Zn"                   "SOM"                  "pH"

## look at the first 6 rows and 7 columns
dat[1:6,1:7]

##   ID Easting Northing Field_No Field_East_Centroid Field_North_Centroid
## 1  1  265632    99300        1            265618.9             99134.54
## 2  2  265625    99275        1            265618.9             99134.54
## 3  3  265650    99275        1            265618.9             99134.54
## 4  4  265600    99250        1            265618.9             99134.54
## 5  5  265625    99250        1            265618.9             99134.54
## 6  6  265650    99250        1            265618.9             99134.54
##   Farmlet_No
## 1          3
## 2          3
## 3          3
## 4          3
## 5          3
## 6          3

## use the sumamry function for the first 7 columns
summary(dat[,1:7])

##        ID            Easting          Northing        Field_No     
##  Min.   :   1.0   Min.   :265425   Min.   :97475   Min.   : 1.000  
##  1st Qu.: 268.2   1st Qu.:265650   1st Qu.:97975   1st Qu.: 3.000  
##  Median : 535.5   Median :265850   Median :98250   Median : 6.000  
##  Mean   : 535.5   Mean   :265907   Mean   :98309   Mean   : 6.313  
##  3rd Qu.: 802.8   3rd Qu.:266175   3rd Qu.:98650   3rd Qu.: 9.000  
##  Max.   :1070.0   Max.   :266528   Max.   :99300   Max.   :15.000  
##  Field_East_Centroid Field_North_Centroid   Farmlet_No  
##  Min.   :265505      Min.   :97629        Min.   :1.00  
##  1st Qu.:265672      1st Qu.:97968        1st Qu.:1.00  
##  Median :265885      Median :98313        Median :2.00  
##  Mean   :265908      Mean   :98308        Mean   :2.01  
##  3rd Qu.:266173      3rd Qu.:98604        3rd Qu.:3.00  
##  Max.   :266451      Max.   :99135        Max.   :3.00

load("data.wk.rda")
## use ls to see what is loaded
ls()

2.2 Saving data

write.csv(dat, file = "new_data.csv")

write.csv(dat, file = "test.csv", row.names = F)

save.image(file = "myworkspace.RData")

# this will save everything in the workspace
save(list = ls(), file = "MyData.RData")
# this will save just dat
save(list = "dat", file = "MyData.RData")
# this will save dat and data.wk
save(list = c("dat", "dtaa.wk"), file = "MyData.RData")

3.1 Initial explorations

We should have a look at the data to see what kinds of properties it has - we could use previously introduced functions such as summary to summarize the attributes. You may wish to use this now. But in terms of model selection what we are interested in the degree to which different variables are good predictors of sulphur.

In order to examine the relationships between the variables in df we can plot the continuous data against each other and use the cor function to examine collinearity amongst them:

plot(df[, -c(1:3)], cex = 0.5, col = grey(0.145,alpha=0.5))

Correlations between numeric variables.

Another useful tool here is to show the upper triangle of the scatterplot matrix with smoothed trend lines. These are achieved with lowess curve fits (Cleveland, 1979) - these are smooth (but possibly curved) bivariate trend lines - and provide a good way of judging by eye whether there is colinearity in the predictor variables. Essentially a straight-line shaped trend with not too much scattering of the individual points suggests collinearity might be an issue. When two predictors are correlated it can be difficult to identify whether it is one or the other (or both) that influence the quantity to be predicted. The code below does this. Note that as collinearity amongst the predictors is the key concern here, we focus on these, ignoring S for this graphic. The upper.panel=panel.smooth causes the lowess curves to be added.

plot(df[,-c(1:4)], cex = 0.5, col = grey(0.145,alpha=0.5), upper.panel=panel.smooth)

Visually, it seems that there may be some colinearity between Ca, pH and P, Slope and P for example.

A numerical table of these relationships as correlations can be generated with the kable function:

library(knitr)
tab <- round(cor(df[,-c(1:3)]), 3)
kable(tab, booktabs = T, format = 'markdown')

What we are looking for in the plots, and to be confirmed in the table, are correlations that look like straight lines, showing that values in one variable change with changes in another.

A final useful tool here is Catherine Hurley’s gclus package. This is similar to the data frame plotting, but provides tools to reorder the variables so that the most strongly correlated pairs are closer to the diagonal in the plots. It can also colour-code the plots to show whether the correlations are positive (red) negative (blue) or near to zero (yellow). For this plot we focus on the predictors again. This is shown below

library(gclus)
cor.mat <- cov(df[,-c(1:3)]) 
cpairs(df[,-c(1:3)],order.single(cor.mat),dmat.color(cor.mat))

The order of the variables in the diagonal generally puts positively correlated variables near to each other as much as possible. The ordering here is based on a clustering method (single linkage) provided via the order.single function in gclus - the colouring is based on correlations (via the dmat.color) function in gclus.

Note, you could also write the table out to a .csv file:

write.csv(cor.mat, file = "tab.csv")

The plot and the table above of the correlations between the variables show that only SOM of the variables are highly correlated (< -0.8 or > +0.8) with each other and with the response variable S and others less so. This provides us with confidence that we can construct a reasonable model of from this data - at least some of the variables are correlated with S - but care may be needed in interpreting some of the regression parameters due to some collinearity in the predictors.

We are interested in a predictive model of income and it is useful to how that value varies. The code below plots the histogram of the variable S. Note the use of the parameter prob=T and prob = F: if set to TRUE this computes probabilities for each histogram bin - rather than a frequency count. This makes the \(y\)-axis compatible with a kernel density estimate (a smooth estimate of the probability distribution of \(x\)) , so the two may be overlaid. Here this is done by using the lines function with a kernel density estimate for the distribution given by density. A rug plot is added to show the exact locations of the values on the \(x\)-axis.

The graphs were obtained with the following code:

hist(df$S, xlab='SulphuR (S)', main='Histogram of Suplhur', col = 'DarkRed')

hist(df$S, prob=T, xlab='SulphuR (S)', main='Histogram of Suplhur', col='NavyBlue')
lines(density(df$S,na.rm=T),col='salmon',lwd=2)
rug(jitter(df$S))

# qqPlot(df$S,main='Normal QQ plot of Sulphur') 

3.2 Regression Models

So we have some data. We have had a look at some of the variables, seen how some of them relate to the variable we are interested in predicting, S, how they correlate with each other, and perhaps we have even started to think about hypotheses that we might wish to test in our model.

We will start with a hypothesized model of Sulphur is associated with Calcium (Ca), Iron (Fe), Phosphorus (P), Soil Organic Matter (SOM), pH, slope (Slope) and landscape aspect (Aspect). In R we can construct a regression model using the function lm which stands for linear model. This is a standard ordinary regression function, although simple it is very powerful and there are many variants and functions for manipulating lm outputs: remember R was originally a stats package.

The equation for a regression is:

\(y_i = \beta_{0} + \sum_{m}^{k=1}\beta_{k}x_{ik} + \epsilon_i \textsf{ for } i \in 1 \cdots n\)

where \(y_i\) is the response variable for observation \(i\), \(x_{ik}\) is the value of the \(k^{\textsf{th}}\) predictor variable for observation \(i\), \(m\) is the number of predictor variables, \(\beta_{0}\) is the intercept term, \(\beta_{k}\) is the regression coefficient for the \(k^{th}\) predictor variable, \(\epsilon_i\) is the random error term for observation \(i\) and \(n\) is the number of observations.

In our case, the response variable \(y\) is sulphur (S), and the 5 chemical variables and 2 landscale variables Ca, Fe, P, SOM, pH, Slope and Aspect are the different predictor variables, the \(k\) different \(x\)’s.

The basic formula in R to specify a linear regression model, lm is as follows:

model <- lm(variable_we_are_trying_to_predict ~ predictor_variable1 +
              predictor_variable2...)

where the ~ sign (tilde) separates the variable we are trying to model (the response variable) from the variables we will use to predict it, the predictor variables. The predictor variables are usually separated by a + sign - meaning that they are linear terms in the model added together.

We can specify our first model.

m <- lm(S~Ca+Fe+P+SOM+pH+Slope+Aspect, data = df) 

And look at its summary attributes

summary(m)

## 
## Call:
## lm(formula = S ~ Ca + Fe + P + SOM + pH + Slope + Aspect, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -333.55  -33.16   -0.40   31.99  285.15 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.466e+02  6.843e+01  10.910  < 2e-16 ***
## Ca           6.213e-02  4.790e-03  12.972  < 2e-16 ***
## Fe          -8.236e-03  3.765e-04 -21.874  < 2e-16 ***
## P            1.152e-01  1.186e-02   9.710  < 2e-16 ***
## SOM          1.952e+01  1.021e+00  19.118  < 2e-16 ***
## pH          -7.060e+01  1.335e+01  -5.289 1.49e-07 ***
## Slope        9.125e+00  8.159e-01  11.184  < 2e-16 ***
## Aspect       1.157e-01  1.481e-02   7.814 1.33e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 51.45 on 1062 degrees of freedom
## Multiple R-squared:  0.8357, Adjusted R-squared:  0.8346 
## F-statistic: 771.5 on 7 and 1062 DF,  p-value: < 2.2e-16

There are lots of things to look at in the summary of the model but we will focus on 3 aspects:

First, Significance. If we examine the model we can see that all are significant predictors of Sulphur. The level of significance is indicated by the \(P\)-value in the column \(Pr(>|t|)\) and the asterisks \(*\). For each coefficient \(\beta_k\) these are the significance levels associated with the hypothesis \(\textsf{H}_0 : \beta_k = 0\) against the alternative \(\textsf{H}_0 : \beta_k \ne 0\) - this corresponds to the probability that one would obtain a \(t\)-statistic for \(\beta_k\) at least as large as the one observed given \(\textsf{H}_0\) is true.

It is a measure of how surprising the relationship is and indicative of whether the relationship between \(x_k\) and \(y\) found in this data could have been found by chance. Very small values suggest that the level of association found here could have come from a random sample of data if \(\beta_k\) really is zero. Cutting to the chase, small \(P\)-values mean strong evidence for a non-zero regression coefficient.

So in this case we can say:

• ‘changes in Calcium are significantly associated with changes in Sulphur at the \(0.001\) level’. This means that the association between the predictor and the response variables is not one that would be found by chance in a series of random samples \(99.999\%\) of the time. Formally we can reject the Null hypothesis that there is no relationship between changes in Sulphur and Calcium
• ‘changes in Fe, P, SOM, pH, Slope and Aspect are significantly associated with changes in S at the at the \(<0.001\) level’, which is actually an indicator of less than \(0.001\). So were are even more confident that the observed relationship between these variables and Sulphur are not by chance and can again confidently reject the Null hypotheses.

Second, the Coeffeicient Estimates. These are the \(\beta_{k}\) values that tell you the rate of change between \(x_k\) and \(y\) in the column ‘Estimate’ in the summary table. The first of these is the Intercept - this is \(\beta_0\) in the regression equation and tells you the value of \(y\) if \(x_1 \cdots x_k\) are all zero.

Remember regression is like taking lots of \(x_k\) and \(y\) values and seeing how well they are correlated. The regression fits a line between them. So in the figure below a regression model with just one \(x\) is defined, it is summarised, and then the \(x\) and \(y\) values are plotted and the model of the relationship is plotted.

m2 <- lm(df$S~df$SOM)
summary(m2)

## 
## Call:
## lm(formula = df$S ~ df$SOM)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -158.12  -50.74  -13.96   36.50  463.37 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 168.2686    10.3170   16.31   <2e-16 ***
## df$SOM       43.1110     0.9616   44.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 74.56 on 1068 degrees of freedom
## Multiple R-squared:  0.653,  Adjusted R-squared:  0.6527 
## F-statistic:  2010 on 1 and 1068 DF,  p-value: < 2.2e-16

plot(df$SOM, df$S, xlab = "x = SOM", ylab = "y = Sulphur")
abline(m2, col = "red", lty = 2, lwd = 2)

Here we can see that the intercept of the relationship is 168.3 and the line shows this, and that the rate of change coefficient between S and SOM, when they are modelled on their own, is NA. This means that each increase in 1 of SOM is associated with an increase in NA in sulphur.

So returning to our model m, we can now see how changes of 1 unit for the different predictor variables are associated with changes in S. The association of P is positive but small, the association of SOM is reduced due to interactions with other \(x\)’s but is still positive, the pH and Fe variables are strongly negatively associated with S.

So the regression model m gives us an idea of how different factors are associated in different ways with sulphur in the study area.

Third the ** \(R^2\) ** value. This provides a measure of model fit. They are calculated as the difference between the actual values of \(y\) (eg of S) and the values predicted by the model - sometimes written as \(\hat{y\ }\) (called ‘\(y\)-hat’). A good description is here http://blog.minitab.com/blog/adventures-in-statistics-2/regression-analysis-how-do-i-interpret-r-squared-and-assess-the-goodness-of-fit: Generally \(R^2\)’s of greater than 0.7 indicate a reasonably well fitting model. Some users of this method describe low \(R^2\) as ‘bad’ - they are not - it just means you a have a more imprecise model - the predictions are a bit wrong-er! Essentially there is a sliding scale of reliability of the predictions - based on the relative amounts of deterministic and random variability in the system being analysed - and \(R^2\) is a measure of this.

Finally, we are also interested in outliers - these are records that the model does not describe very well, either vastly over-predicting Sulphur or under-predicting it. When the bivariate regression was considered earlier, these were easy to find, since they were just the points lying a long way from the regression line. However, when there are multiple predictors, relationships cannot be encapsulated in two dimensional space.

For this situation, a useful diagnostic is the studentised residual. For different values of predictors, the expected difference between the actual and fitted \(y\) values has a different variance even if the theoretic model gives \(\epsilon_i\) a constant variance. This discrepancy is due to issues such as collinear predictors. Standardised variables attempt to compensate for this. Plotting these and looking for notably high or low values is a helpful way of spotting outliers in the regression model.

A simple way of doing this is just to plot the residuals in order:

s.resids <- rstudent(m)
plot(s.resids,type='h')
abline(h=c(-2,2),col='lightskyblue')

Here we see some residuals are relatively high. Normally, one might expect standardised residuals to lie between -2 and 2. Adding lines to this plot helps to identify when this is not the case - as is done above.

Another useful job is to identify which are the outlying cases - here just done as a basic table (via kable)

df_out <- data_frame(Name=df$ID,Residual=s.resids)
df_out <- df_out[abs(df_out$Residual)>2,]
kable(df_out)

Negative values here have a much lower income than the model predicts, positive values suggest it is much higher.

Finally, we need to check another model assumption in lm - that the error term \(\epsilon_i\)’s are normally distributed. The best way to do this is via a qq-plot which plots the empirical quantiles of a variable whose distribution we are interested in against the theoretical distribution we are testing for. If the variable is from this distribution, the results should yield a straight line - frequently a reference line is provided to test this. Here s.resids is the variable of interest - the studentised residuals, and these are compared to the \(t\) distribution. Although the \(\epsilon_i\)’s have a normal distribution, theory shows that the studentised residuals derived from a sample coming from this model have a \(t\) distribution. However, the qqPlot function from the package car does all the work, and requires only the regression model object to create a plot:

qqPlot(m)

This demonstrates that mostly the residuals do conform to this, with a couple of extreme outliers at the higher end - labelled in the plot. Consulting the earlier table suggests these are Fayette and ID 1009, 1047, 465, and 1070

It is also important to consider spatial patterns in the residuals, from a geographical viewpoint. The code below plots the outliers and their spatial distributions:

# determine residuals
df$s.resids = rstudent(m)

# 1.5 standard deviations
resid.shades = shading(c(-1.5,1.5),c("indianred1","antiquewhite2","powderblue"))
# define a color scheme
cols = resid.shades$cols[1 + findInterval(df$s.resids, resid.shades$breaks)]
# set plot paramters
par(mar = c(0,0,2,0))

# make some spatial data 
# define a projection
os.proj <- CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs ")

df.sp <- SpatialPointsDataFrame(coords = df[, 2:3], data = data.frame(df), 
                                proj4string = os.proj)

# have a look at the spatial distribution of outliers
choropleth(df.sp,df$s.resids,resid.shades, main = "Model Outliers", pch = 19, cex = 0.7) 

The regression session finishes here, but you will need the results of this session later on. To save all of the variables, regression objects and so on from this session, enter

save.image(file='session2.RData')

This saves an image of all of the R objects currently defined in a file called session2.RData - later on you will be able to load this again via the load function.

4.1 Reading and Writing Spatial Data

As we are working with spatial data it seems appropriate to consider how to get spatial data in and out of R.

Very often we have data that is in a particular format such as shapefile format. R has the ability to read and write data from and to many different spatial data formats using functions in the rgdal and sf packages - we will consider them both here:

However, as the sp package and data formats are depreciated, so are their functions for reading and writing spatial data. The rgdal package includes 2 generic functions for reading and writing all kinds of spatial data - pretty much you name the format and it has a driver for that format. These are

readOGR()
writeOGR()

Spatial data can be read into R and assigned to a variable using the readOGR function:

fields <- readOGR("FP_Fields_Fenced_post13082013.shp") 

## OGR data source with driver: ESRI Shapefile 
## Source: "FP_Fields_Fenced_post13082013.shp", layer: "FP_Fields_Fenced_post13082013"
## with 20 features
## It has 13 fields

If you enter:

class(fields)

You will see that the class of the fields object is sp.

The fields sp format object - we can write this out to a shape file using the writeOGR() function as follows:

writeOGR(obj=fields, dsn=".", layer="fields", driver="ESRI Shapefile") 

You will see that a shapefile has been written into your current working directory, with its associated supporting files (.dbf, etc) that can be recognised by other applications (QGIS etc). You should examine the writeOGR and readOGR function in the rgdal package. R is also able to read and write other proprietary spatial data formats using a number of packages, which you should be able to find through a search of the R help system or via an internet search engine. The rgdal package, which is the R version of the Geospatial Data Abstraction Library. It includes a number of methods for reading and writing spatial objects, including to and from SpatialXDataFrame objects:

The full syntax can be important - the code below overwrites any existing similarly named file:

writeOGR(fields, dsn = ".", layer = "fields", 
         driver="ESRI Shapefile", overwrite_layer = T)

The dsn parameter is important here: for shape files it determines the folder the files are written to. In the above example it was set to "." which places the files in the current working directory. You could specify a file path here.

For a PC might be something like:D:\MyDocuments\MyProject\DataFiles
And for a MAC: /Users/lex/my_docs/project

The setwd() and getwd() functions can be used in determining and setting the file path. You may want to set the file path and then use the dsn setting as above:

setwd("/Users/lex/my_docs/project")
writeOGR(fields, dsn = ".", layer = "fields", 
         driver="ESRI Shapefile", overwrite_layer = T)

Or you could use the getwd() function, save the results to variable and pass this to writeOGR:

td <- getwd()
writeOGR(fields, dsn = td, layer = "fields", 
         driver="ESRI Shapefile", overwrite_layer = T)

You should also examine the functions for reading and writing raster layers in rgdalwhich are readGDAL and writeGDAL. These read and write functions in rgdal are incredibly powerful and can read / write almost any format of spatial data.

Spatial data can be also be read in and written out using the sf functions:

st_read()
st_write()

They work in a similar way - we do not have the time to examine them in detail here.

4.2 Mapping with tmap

The tmap package is for thematic visualization, with a focus on mapping the spatial distribution of data attributes. It has a similar grammar to plotting with ggplot in that it seeks to handle each element of the map separately in a series of layers and in so doing to exercise control over each element. This is different to the basic plot functions used above to map sp and sf data.In earlier sections qtm() was used to compose a quick map and now we will expand this more fully. The code in this section will use the fields and df.sp datasets.

A number of plot parameters additional to the ones that have previously been used exist, including different window sizes, multiple plots in the same window, polygon or area shading, hatching, boundary thickness, boundary colour, labeling etc.

A useful introductory overview of tmap can be found in the tmap-nutshell vignette. You can uses this later to get some further ideas about mapping options:

vignette("tmap-nutshell")

Further aesthetics for derived layers are described in

help(tmap)

The code below recreates and augments the map above created with the choropleth function in GISTools this time with filed boundaries.

tm_shape(df.sp) +
  tm_dots("s.resids", size = 0.1, breaks = c(-10,-1.5,2,10), palette = "RdBu") +
  tm_shape(fields) + 
  tm_borders()

So what you can see in the above code are two sets of tmap plot commands: the first set plots the df.sp dataset, specifying a dots to show the data locations, a size and some breaks. The second set adds the field outline.

It is also possible to plot multiple different maps from different datasets together but this requires a bit more control over the tmap parameters.

One of the great features of tmap is its ability to include leaflet plots. This can be done by setting the tmap_mode parameter:

tmap_mode('view')

## tmap mode set to interactive viewing

tm_shape(df.sp) +
  tm_dots("s.resids", breaks = c(-10,-1.5,2,10), palette = "RdBu") +
  tm_shape(fields) + 
  tm_borders()

Try changing the layers icon - the Esri.WorldTopoMap options! This is where Harry works - on a little farm at the end of the world!!! Note the street name of Cocktree Throat - I thought this was an illness not a place!

Mapping GW outputs

The variation in the GWR coefficient estimates can be mapped to show how different factors are varying associated with Sulphur in different parts of the study area

Maps where large variation in the relationship between the predictor variable and sulphur was found:

tmap_mode('plot')

## tmap mode set to plotting

tm_shape(m4$SDF) +
  tm_dots(c("SOM", "Slope"), size = 0.1) +
  tm_layout(legend.position = c("right","top"))

Maps showing coefficients who values flip from positive to negative in different parts of the study area

tm_shape(m4$SDF) +
  tm_dots(c("Slope", "pH"), size = 0.1) +
  tm_layout(legend.position = c("right","top")) +
  tm_style_col_blind()

So what we see here are the local variations in the degree to which changes in different variables are associated with changes in Sulphur The maps above are not simply maps of the predictor variable.