image.png

<aside> 📦

You have to first install and load these packages in R

</aside>

Handling vector data with sf


  1. Our first exercise is to draw the famous and beautiful Wytham Woods near us, which hosts the world’s longest ecological experiment. Wytham Woods is a semi-natural woodland, where different parts are established at different times. The stands can be classified as N (natural), M (mixed), P-E19 (early 19th-century plantation), P-M19 (mid 19th-century plantation), and P-M20 (mid 20th-century plantation).

    <aside> 📥

    Download

    ww.zip

    oaks.csv

    </aside>

    1. Use st_read to load the correct shapefile, save it as ww.
    2. Call ww. How many features does it store?
    3. Similar to a standard dataframe, we can subset ww (a sf object) by using the square bracket notation and use the drop argument to drop geometries.
      1. ww[1,]
      2. ww[nc$stand == "P-M19", ]
      3. ww[1, nc$stand == "P-M19", ]
      4. ww[1, nc$stand == "P-M19", drop = TRUE]
      5. [st_geometry](<https://r-spatial.github.io/sf/reference/st_geometry.html>)(ww)[[1]]
    4. We can use plot(ww) to quickly visualise the map of the first attributes.
      1. We can also specify only the strand attribute by the square bracket notation.
    5. We will use ggplot2 to produce a fancier, publication-ready figure.
      1. Remember that ggplot2 uses layer-based aesthetics, we will start with an empty grid ggplot(ww).
      2. We need to add the sf object onto it, and specify to use stand to colour code. Hint: use geom_sf() and aes().
      3. We will then add a layer called coord_sf() to make sure the coordinate grid can be added.
      4. Complete your figure with a title and a legend, and perhaps you can also use manual colour palette that might match the stand type.
      • Answer
    6. We will produce a choropleth by calculating the area for each stand.
      1. Use mutate and st_read to calculate the area for each stand, name it as area_m2.
      2. By default, the units (m^2) will also be stored in the variable, but they will not be understood by other functions such as ggplot. It might sound counterproductive, but use drop_units to get rid of them.
      3. Draw a map with different colours corresponding to area of each stand.
      • Answer
    7. Now, investigate the oaks.csv. It is a simple CSV with latitudinal and longitudinal coordinates.
      1. Read in the data frame and convert it to sf points, using st_as_sf. Remember to set the correct corrds = and crs = .
      2. Similar to the ggplot logic above, use geom_sf to add the trees onto the map of Wytham Woods as black dots.
      3. Use st_intersect to count the number of points with the polygons. What is the returned object?
      4. Then use the lengths to calculate the number of points.
      5. Now, draw another map that shows the number of trees per stand.
      • Answer

Handling raster and vector data with terra


  1. We will explore the temperature data across UK with the famous and useful WorldClim database.

    <aside> 📥

    Download

    tavg.tif

    prec.tif

    srad.tif

    wind.tif

    </aside>

    <aside> 💡

    You can actually use the geodata package to directly download WordClim data like this:

    geodata::worldclim_country(country = "United Kingdom", var = "tavg", res = 2.5, path = "."))
    

    </aside>

    1. Use rast() to read the raster and save it as tavg.

    2. Inspect the objects:

      1. What is the class of tavg and ww_v?
      2. Use res(), ext(), and crs() on tavg. What do they tell you?
      3. Use summary(ww.v) to inspect the attribute table.
    3. Use plot() to quickly visualise the temperature raster.

    4. It is a Raster Stack which saves monthly variable. Let’s use mean(tavg) to create an annual average tavg.ann. Visualise it again.

    5. We can then also download the map of UK with the rnaturalearth package (South 2017), with the function ne_states.

      uk <- rnaturalearth::ne_states("United Kingdom", returnclass = "sf")
      
    6. Use ggplot() + geom_sf() to visualise this UK map.

    7. You will actually notice that the extent of tavg.ann raster and the uk map is slightly different. We will use crop() to crop the raster to the extent of uk.

      • Answer
    8. There’s a bit of France in tavg.ann. We can use mask() to convert all values outside the uk map to NA.

      • Answer
    9. The aggregate() function can be used to aggregate groups of cells of a raster in order to create a new raster with a lower resolution (thus larger cells). It takes two arguments. fact = is the aggregation factor expressed as number of cells in each direction (both horizontally and vertically), or two integers denoting the horizontal and vertical aggregation factor separately. fun = specifies the function used to aggregate values.

      1. Aggregate the mean by an overall factor of 10.
      2. Aggregate the mean by an overall factor of 50.
      3. Aggregate the mean by an vertical factor of 50, and no change in horizontal direction.
      4. Aggregate the maximum value by an overall factor of 10.
      • Answer
    10. We can use centroids() to calculate the centroid of each district in the UK (each polygon), but we have to convert the uk to a vectorised form with vect(). It is now a SpatVector object**.**

    11. We can then extract the raster values at a set of points using the extract() function. Here we will extract the mean annual temperature of the centroid of each district.

      1. First, plot tavg.ann again.
      2. Second, plot vect(uk) with add = T.
      3. Third, use crds() to extract the coordinates of those centroids.
      4. Fourth, add them as points().
      5. Now we can use extract() to pass the SpatVector object with the raster as the first argument, then a data frame of centroids as the second argument.
      • Answer
    12. We can also use the same extract() function to obtain the average raster values of tavg.ann by polygons in uk. We will save this as a feature inside uk.

      uk$tavg.ann <- extract(tavg.ann, vect(uk), mean, na.rm = T)$mean
      ggplot(uk) + geom_sf(aes(fill = tavg.ann))
      
    13. We could also have computed the average by weighting the area of each polygon, simply set weights = T.

      uk$tavg.ann.weighted <- extract(tavg.ann, vect(uk), mean, weights = T, na.rm = T)$mean
      ggplot(uk) + geom_sf(aes(fill = tavg.ann.weighted))
      

<aside> 🎉

Congratulations! You now have learned how to handle basic manipulation and visualisation of vectors and rasters with the two most important packages sf and terra. Drawing a map in R is no longer a challenge!

</aside>