Chapter 3 Rasters, descriptive statistics and interpolation

3.1 Learning outcomes

By the end of this practical you should be able to:

  1. Load, manipulate and interpret raster layers
  2. Observe and critique different descriptive data manipulation methods and outputs
  3. Execute interpolation of points to a raster layer
  4. Construct a methodology for comparing raster

3.2 Homework

Outside of our schedulded sessions you should be doing around 12 hours of extra study per week. Feel free to follow your own GIS interests, but good places to start include the following:

Assignment

This week you need to assess the validity of your topic in terms of data suitability, then sketch out an introduction and literature review including academic and policy documents (or any reputable source) — it’s the same as last week!

Reading

This week:

Remember this is just a starting point, explore the reading list, practical and lecture for more ideas.

3.4 Introduction

This practical is composed of three parts. To start with we’re going to load some global raster data into R. In the second part we extract data points (cities and towns) from this data and generate some descriptive statistics and histograms. In the final section we explore interpolation using point data.

3.5 Part 1 rasters

So far we’ve only really considered vector data. Within this practical we will explore some raster data sources and processing techniques. If you recall rasters are grids of cell with individual values. There are many, many possible sources to obtain raster data from as it is the data type used for the majority (basically all) of remote sensing data.

3.5.1 WorldClim data

To start with we are going to use WorldClim data — this is a dataset of free global climate layers (rasters) with a spatial resolution of between 1\(km^2\) and 240\(km^2\).

  1. Download the data from: http://worldclim.org/version2

  2. Select any variable you want at the 5 minute second resolution.

What is a 5 minute resolution i hear you ask? Well, this geographic reference system treats the globe as if it was a sphere divided into 360 equal parts called degrees. Each degree has 60 minutes and each minute has 60 seconds. Arc-seconds of latitude (horizontal lines in the globe figure below) remain almost constant whilst arc-seconds of longitude (vertical lines in the globe figure below) decrease in a trigonometric cosine-based fashion as you move towards the Earth’s poles. This causes problems as you increase or decrease latitude the longitudial lengths alter…For example at the equator (0°, such as Quito) a degree is 111.3 km whereas at 60° (such as Saint Petersburg) a degree is 55.80 km …In contrast a projected coordinate system is defined on a flat, two-dimensional plane (through projecting a spheriod onto a 2D surface) giving it constant lengths, angles and areas…

This figure is taken directly from Lovelace et al. (2019) section 2.2. Illustration of vector (point) data in which location of London (the red X) is represented with reference to an origin (the blue circle). The left plot represents a geographic CRS with an origin at 0° longitude and latitude. The right plot represents a projected CRS with an origin located in the sea west of the South West Peninsula.This figure is taken directly from Lovelace et al. (2019) section 2.2. Illustration of vector (point) data in which location of London (the red X) is represented with reference to an origin (the blue circle). The left plot represents a geographic CRS with an origin at 0° longitude and latitude. The right plot represents a projected CRS with an origin located in the sea west of the South West Peninsula.

Figure 3.1: This figure is taken directly from Lovelace et al. (2019) section 2.2. Illustration of vector (point) data in which location of London (the red X) is represented with reference to an origin (the blue circle). The left plot represents a geographic CRS with an origin at 0° longitude and latitude. The right plot represents a projected CRS with an origin located in the sea west of the South West Peninsula.

If you are still a bit confused by coordiate reference systems then stop and take some time to have a look at the resources listed here. It is very important to understand projection systems.

This is the best resources i’ve come across explaining coordiate reference systems are:

Others include:

This YouTube video that shows the differences between geographic and projected coordinate reference systems and some limitations of ArcMap…

And this great YouTube produced by SciShow shows the comporsises of different map projections

  1. Unzip and move the data to your project folder. Now load the data. We could do this individually….
## class      : RasterLayer 
## dimensions : 2160, 4320, 9331200  (nrow, ncol, ncell)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/Andy/OneDrive - University College London/Teaching/CASA0005repo/prac3_data/wc2.0_5m_tavg_01.tif 
## names      : wc2.0_5m_tavg_01 
## values     : -46.697, 34.291  (min, max)
  1. Then have a quick look at the data

Now we can actually see some data…here is a quick example of using the Robinson projection saved to a new variable. Don’t worry about the code, just take in that it is different to the above plot

  1. A better and more efficient way is to firstly list all the files stored within our directory
##  [1] "prac3_data/wc2.0_5m_tavg_01.tif" "prac3_data/wc2.0_5m_tavg_02.tif"
##  [3] "prac3_data/wc2.0_5m_tavg_03.tif" "prac3_data/wc2.0_5m_tavg_04.tif"
##  [5] "prac3_data/wc2.0_5m_tavg_05.tif" "prac3_data/wc2.0_5m_tavg_06.tif"
##  [7] "prac3_data/wc2.0_5m_tavg_07.tif" "prac3_data/wc2.0_5m_tavg_08.tif"
##  [9] "prac3_data/wc2.0_5m_tavg_09.tif" "prac3_data/wc2.0_5m_tavg_10.tif"
## [11] "prac3_data/wc2.0_5m_tavg_11.tif" "prac3_data/wc2.0_5m_tavg_12.tif"
  1. Then load all of the data straight into a raster stack. A raster stack is a collection of raster layers with the same spatial extent and resolution.
## class      : RasterStack 
## dimensions : 2160, 4320, 9331200, 12  (nrow, ncol, ncell, nlayers)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : wc2.0_5m_tavg_01, wc2.0_5m_tavg_02, wc2.0_5m_tavg_03, wc2.0_5m_tavg_04, wc2.0_5m_tavg_05, wc2.0_5m_tavg_06, wc2.0_5m_tavg_07, wc2.0_5m_tavg_08, wc2.0_5m_tavg_09, wc2.0_5m_tavg_10, wc2.0_5m_tavg_11, wc2.0_5m_tavg_12 
## min values :          -46.697,          -44.559,          -57.107,          -62.996,          -63.541,          -63.096,          -66.785,          -64.600,          -62.600,          -54.400,          -42.000,          -45.340 
## max values :           34.291,           33.174,           33.904,           34.629,           36.312,           38.400,           43.036,           41.073,           36.389,           33.869,           33.518,           33.667

In the raster stack you’ll notice that under dimensions there are 12 layers (nlayers). The stack has loaded the 12 months of average temperature data for us in order.

  1. To access single layers within the stack:
## class      : RasterLayer 
## dimensions : 2160, 4320, 9331200  (nrow, ncol, ncell)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/Andy/OneDrive - University College London/Teaching/CASA0005repo/prac3_data/wc2.0_5m_tavg_01.tif 
## names      : wc2.0_5m_tavg_01 
## values     : -46.697, 34.291  (min, max)
  1. We can also rename our layers within the stack:
  1. Now to get data for just January use our new layer name
## class      : RasterLayer 
## dimensions : 2160, 4320, 9331200  (nrow, ncol, ncell)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/Andy/OneDrive - University College London/Teaching/CASA0005repo/prac3_data/wc2.0_5m_tavg_01.tif 
## names      : Jan 
## values     : -46.697, 34.291  (min, max)

3.6 Part 2 descriptive statistics

Descriptive (or summary) statistics provide a summary of our data, often forming the base of quantitiatve analysis leading to inferential statistics which we use to make infereces about our data (e.g. judegements of the probability that the observed difference between two datasets is not by chance)

3.6.1 Data preparation

  1. Let’s take Perth as an example. We can subset our data either using the row name:
  1. Or the row location:

3.6.2 Histogram

A histogram lets us see the frequency of distribution of our data.

  1. Make a histogram of Perth’s temperature

Remember what we’re looking at here. The x axis is the temperature and the y is the frequency of occurrence.

  1. That’s a pretty simple histogram, let’s improve the aesthetics a bit.

  1. Check out the histogram information R generated
## $breaks
## [1] 12 14 16 18 20 22 24 26
## 
## $counts
## [1] 1 3 2 2 1 1 2
## 
## $density
## [1] 0.04166667 0.12500000 0.08333333 0.08333333 0.04166667 0.04166667 0.08333333
## 
## $mids
## [1] 13 15 17 19 21 23 25
## 
## $xname
## [1] "Perthtemp"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

Here we have:

  • breaks — the cut off points for the bins (or bars), we just specified these
  • counts — the number of cells in each bin
  • midpoints — the middle value for each bin
  • density — the density of data per bin

3.6.3 Using more data

This was still a rather basic histogram, what if we wanted to see the distribution of temperatures for the whole of Australia in Jan (from averaged WorldClim data) as opposed to just our point for Perth.

  1. First, we need to source and load a vector of Australia. Go to: https://gadm.org/download_country_v3.html and download the GeoPackage

  2. Check what layers are within a GeoPackage using:

## Driver: GPKG 
## Available layers:
##     layer_name geometry_type features fields
## 1 gadm36_AUS_0 Multi Polygon        1      2
## 2 gadm36_AUS_1 Multi Polygon       11     10
## 3 gadm36_AUS_2 Multi Polygon      569     13
  1. Then read in the GeoPackage layer for the whole of Australia
## Reading layer `gadm36_AUS_0' from data source `C:\Users\Andy\OneDrive - University College London\Teaching\CASA0005repo\prac3_data\gadm36_AUS.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 112.9211 ymin: -55.11694 xmax: 159.1092 ymax: -9.142176
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Check the layer by plotting the geometry…we could do this through…

But as the .shp is quite complex (i.e. lots of points) we can simplify it first with the rmapshaper package — install that now..if it doesn’t load (or crashes your PC) this isn’t an issue. It’s just good practice that when you load data into R you check to see what it looks like…

## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 112.9211 ymin: -54.7775 xmax: 159.1054 ymax: -9.221099
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   GID_0    NAME_0                           geom
## 1   AUS Australia MULTIPOLYGON (((158.8517 -5...

This should load quicker, but for ‘publication’ or ‘best’ analysis (i.e. not just demonstrating or testing) i’d recommend using the real file to ensure you don’t simply a potentially important variable.

Check out this vignette for more information about rmapshaper

  1. Next, set our map extent to the outline of Australia then crop our WorldClim dataset to it
## class      : Extent 
## xmin       : 112.9211 
## xmax       : 159.1092 
## ymin       : -55.11694 
## ymax       : -9.142176
## class      : RasterBrick 
## dimensions : 551, 554, 305254, 12  (nrow, ncol, ncell, nlayers)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : 112.9167, 159.0833, -55.08333, -9.166667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      :       Jan,       Feb,       Mar,       Apr,       May,       Jun,       Jul,       Aug,       Sep,       Oct,       Nov,       Dec 
## min values :  6.360000,  6.287500,  5.554286,  4.051429,  2.542857, -0.171000, -1.945000, -1.677000,  0.677000,  3.054286,  3.842424,  5.433333 
## max values :  34.29100,  33.17400,  32.39700,  30.07200,  28.50000,  27.40000,  26.90000,  27.20000,  29.31209,  31.72000,  33.51800,  33.66700

You’ll notice that whilst we have the whole of Australia the raster hasn’t been perfectly clipped to the exact outline….the extent just specifies an extent box that will cover the whole of the shape.

  1. If want to just get raster data within the outline of the shape:

You could also run this using the original worldclimtemp raster, however, it may take some time. I’d recommend cropping to the extent first.

Both our Austemp and exactAus are raster bricks. A brick is similar to a stack except it is now stored as one file instead of a collection.

  1. Let’s re-compute our histogram for Australia in March. We could just use hist like we have done before

However we have a bit more control with ggplot()

3.6.4 Histogram with ggplot

  1. We need to make our raster into a data.frame to be compatible with ggplot2

How about plotting multiple months of temperature data on the same histogram

  1. As we did in practical 2, we need to put our variaible (months) into a one coloumn using melt. We will do this based on the names of our coloumns in alldf…

We could also use the functions we learnt about in Tidying data

  1. Then subset the data, selecting two months
  1. Get the mean for each month we selected
##   Month grp.mean
## 1   Jan 28.11321
## 2   Jun 14.96415
  1. Select the colour and fill based on the variable (which is our month). The intercept is the mean we just calculated, with the lines also based on the coloumn variable.

Note how i adjusted the title after i selected the theme, if i had done this before the theme defaults would have overwritten my command.

  1. Have you been getting an annoying error message about bin size and non-finate values? Me too!…Bin size defaults to 30 in ggplot2 and the non-finate values is referring to lots of NAs (no data) that we have in our dataset. In the code below i’ve selected a bin width of 5 and removed all the NAs with complete.cases and produced a faceted plot…
## [1] 1201812       2
## [1] 3663048       2

Does this seem right to you? Well…yes. It shows that the distribution of temperature is higher (or warmer) in the Australian summer (Dec-Feb) than the rest of the year, which makes perfect sense.

How about an interactive histogram using plotly

  1. See if you can understand what is going on in the code below. Run each line seperately.

This format of code where you set lots of varaibles then call them within a plot, package or fuction is sometihng you should become more familiar with as it’s considerd good practice. If you were to go on and produce multiple plots using the same legends / aesthetics you only ahve to set them once.

Ok so enough with the histograms…the point is to think about how to best display your data both effectively and efficiently.

  1. Let’s change the pace a bit and do a quickfire of other descrptive statistics you might want to use…
##   variable grp.mean
## 1      Jan 28.11321

3.7 Part 3 interpolation

What if you had a selection of points over a spatial area but wanted to generate a complete raster. For this example, we will take our sample points (Australian cities) and estimate data between them using interpolation.

  1. If you look at our samples and AUcitytemp data the lat and lon is only in the former. We need to have this with our temperature data so let’s combine it using cbind
  1. Now we need to tell R that our points are spatial points
  1. You’ll notice that here i’ve just nicked the CRS from our worldclimtemp. In generally it’s good practice to avoid using static or hard coding references, by that i mean if we added another coloumn to our samplestemp data (or manipulated somehow) then using this…don’t run this….

This would give us a headache as coloumns 13 and 14 might no longer be longitude and latitude.

  1. Right…plot the Australian geometry outline then add our spatial data points ontop…

  1. Let’s interpolate using Inverse Distance Weighting, or IDW as it’s more commonly known. IDW is a deterministic method for multivaraite interpolation that estaimtes values for a surface using a weighted average of the provided data. The values closer to the point being predicted have more weight than those further away. The rate at which distance from the provided point imapcts the predcted point is controlled by the power of p. If p=0 then there is no decrease with distance.

For more infomation see: https://pro.arcgis.com/en/pro-app/help/analysis/geostatistical-analyst/how-inverse-distance-weighted-interpolation-works.htm

  1. To get a meaningful result we could run some more calucaltions on let’s project our data to GDA94 (EPSG:3112)

I’ve been a bit lazy here, if you recall from the previous practical session there are different ways to reproject data based if you have an SP or SF object. SP objects require you to specify the entire CRS string, so i just convered it to SF, used the ESPG code and converted it back to SP. The main reason for doing this is that i made a grid to store my interpolation and having a remote sensing background i wanted to specify the pixel size. The equivalent function if SF won’t let you specify pixel size or there is no easy and straightforward way to do it (at least to my knowledge).

  1. Next, create an empty grid where cellsize is the spatial resolution, cellsize will overwrite the number of pixels we specified (n). Here as we’ve used a projected CRS i’ve put a high cellsize (in metres) so 200km by 200km cells. You can use a smaller number if you wish but it will take much longer to process.
## [inverse distance weighted interpolation]

IDW is just one method for interpolating data, there are many more, if you are interested check out: https://mgimond.github.io/Spatial/interpolation-in-r.html

3.8 Auto data download

In this practical I’ve shown you how to source the data online, download it and load it into R. However for both WorldClim and GADM we can do this straight from R using the getData function….i’m sorry for making you do it the long way, but it’s good to do things manually to see how they work.

WARNING, this may take some time. I’ve changed the resolution to 10 degrees, but I’d advise not running this in the practical session.

Now for GADM

Much more convenient right?

3.9 Advanced analysis

Are you already comptent with raster analysis and R, then have a go at completing this task in the practical session.

Within the practical we’ve loaded one and created one raster layer. Undertake some comparative analysis to detemrine spatial (and temporal if appropraite) differences between the rasters here and any others you may wish to create (e.g. from other interpolation methods). Try to identify where the varaitions are and explain why they are occuring.

You could assume that one raster is the ‘gold standard’ meaning it’s beleived to be fully correct and compare others to it.

… Or you could go further than this and obtain weather station temperature data (or any other variable) for multiple sites, interpolate based on 50% of the sites and use the remaining sites to assess the accuracy of your selected method / the WorldClim data.

Free weather station data can be found here: https://rp5.ru/Weather_in_the_world

Have a go and discuss with your fellow students / members of the teaching team during the practical sessions or on slack.

3.10 Feedback

Was anything that we explained unclear this week or was something really clear…let us know using the feedback form. It’s anonymous and we’ll use the responses to clear any issues up in the future / adapt the material.