diff --git a/01-intro-to-r.md b/01-intro-to-r.md index 1b99653e..ebc3ec6c 100644 --- a/01-intro-to-r.md +++ b/01-intro-to-r.md @@ -84,7 +84,7 @@ An alternative solution is to create the folders using R command `dir.create()`. In the console type: -```r +``` r dir.create("data") dir.create("data_output") dir.create("documents") @@ -197,7 +197,7 @@ We will however need to install the `here` package. To do so, please go to your script and type: -```r +``` r install.packages("here") ``` @@ -210,7 +210,7 @@ of installed packages. If not, you will need to install it by writing in the script: -```r +``` r install.packages('tidyverse') ``` @@ -231,7 +231,7 @@ execution. Thanks to this feature, you can annotate your code. Let's adapt our script by changing the first lines into comments: -```r +``` r # install.packages('here') # install.packages('tidyverse') ``` @@ -243,7 +243,7 @@ Installing packages is not sufficient to work with them. You will need to load them each time you want to use them. To do that you use `library()` command: -```r +``` r # Load packages library(tidyverse) library(here) @@ -275,7 +275,7 @@ It might be confusing, so let's see how it works. We will use the `here()` funct from the `here` package. In the console, we write: -```r +``` r here() here('data') ``` @@ -292,7 +292,7 @@ We will save it in the `data/` folder, where the **raw** data should go. In the script, we will write: -```r +``` r # Download the data download.file( "https://bit.ly/geospatial_data", @@ -321,7 +321,7 @@ will not cover these in the workshop. You can use R as calculator, you can for example write: -```r +``` r 1 + 100 1 * 100 1 / 100 @@ -335,14 +335,14 @@ use them whenever we need to. We using the assignment operator `<-`, like this: -```r +``` r x <- 1 / 40 ``` Notice that assignment does not print a value. Instead, we've stored it for later in something called a variable. `x` variable now contains the value `0.025`: -```r +``` r x ``` @@ -352,21 +352,21 @@ Our variable `x` can be used in place of a number in any calculation that expect a number, e.g. when calculating a square root: -```r +``` r sqrt(x) ``` Variables can be also reassigned. This means that we can assign a new value to variable `x`: -```r +``` r x <- 100 x ``` You can use one variable to create a new one: -```r +``` r y <- sqrt(x) # you can use value stored in object x to create y y ``` diff --git a/02-data-structures.md b/02-data-structures.md index 950480a5..319a9995 100644 --- a/02-data-structures.md +++ b/02-data-structures.md @@ -64,33 +64,33 @@ Note that vector data in the geospatial context is different from vector data ty You can create a vector with a `c()` function. -```r +``` r # vector of numbers - numeric data type. numeric_vector <- c(2, 6, 3) numeric_vector ``` -```output +``` output [1] 2 6 3 ``` -```r +``` r # vector of words - or strings of characters- character data type character_vector <- c('Amsterdam', 'London', 'Delft') character_vector ``` -```output +``` output [1] "Amsterdam" "London" "Delft" ``` -```r +``` r # vector of logical values (is something true or false?)- logical data type. logical_vector <- c(TRUE, FALSE, TRUE) logical_vector ``` -```output +``` output [1] TRUE FALSE TRUE ``` @@ -99,21 +99,21 @@ logical_vector The combine function, `c()`, will also append things to an existing vector: -```r +``` r ab_vector <- c('a', 'b') ab_vector ``` -```output +``` output [1] "a" "b" ``` -```r +``` r abcd_vector <- c(ab_vector, 'c', 'd') abcd_vector ``` -```output +``` output [1] "a" "b" "c" "d" ``` @@ -141,27 +141,27 @@ A common operation you want to perform is to remove all the missing values (in R denoted as `NA`). Let's have a look how to do it: -```r +``` r with_na <- c(1, 2, 1, 1, NA, 3, NA ) # vector including missing value ``` First, let's try to calculate mean for the values in this vector -```r +``` r mean(with_na) # mean() function cannot interpret the missing values ``` -```output +``` output [1] NA ``` -```r +``` r # You can add the argument na.rm=TRUE to calculate the result while # ignoring the missing values. mean(with_na, na.rm = T) ``` -```output +``` output [1] 1.6 ``` @@ -171,22 +171,22 @@ For this you need to identify which elements of the vector hold missing values with `is.na()` function. -```r +``` r is.na(with_na) # This will produce a vector of logical values, ``` -```output +``` output [1] FALSE FALSE FALSE FALSE TRUE FALSE TRUE ``` -```r +``` r # stating if a statement 'This element of the vector is a missing value' # is true or not !is.na(with_na) # The ! operator means negation, i.e. not is.na(with_na) ``` -```output +``` output [1] TRUE TRUE TRUE TRUE FALSE TRUE FALSE ``` @@ -195,14 +195,14 @@ Now we need to retrieve the subset of the `with_na` vector that is not `NA`. Sub-setting in `R` is done with square brackets`[ ]`. -```r +``` r without_na <- with_na[ !is.na(with_na) ] # this notation will return only # the elements that have TRUE on their respective positions without_na ``` -```output +``` output [1] 1 2 1 1 3 ``` @@ -224,22 +224,22 @@ Once created, factors can only contain a pre-defined set of values, known as levels. -```r +``` r nordic_str <- c('Norway', 'Sweden', 'Norway', 'Denmark', 'Sweden') nordic_str # regular character vectors printed out ``` -```output +``` output [1] "Norway" "Sweden" "Norway" "Denmark" "Sweden" ``` -```r +``` r # factor() function converts a vector to factor data type nordic_cat <- factor(nordic_str) nordic_cat # With factors, R prints out additional information - 'Levels' ``` -```output +``` output [1] Norway Sweden Norway Denmark Sweden Levels: Denmark Norway Sweden ``` @@ -251,19 +251,19 @@ This can come in handy when performing statistical analysis. You can inspect and adapt levels of the factor. -```r +``` r levels(nordic_cat) # returns all levels of a factor vector. ``` -```output +``` output [1] "Denmark" "Norway" "Sweden" ``` -```r +``` r nlevels(nordic_cat) # returns number of levels in a vector ``` -```output +``` output [1] 3 ``` @@ -281,7 +281,7 @@ displayed in a plot or which category is taken as a baseline in a statistical mo You can reorder the categories using `factor()` function. This can be useful, for instance, to select a reference category (first level) in a regression model or for ordering legend items in a plot, rather than using the default category systematically (i.e. based on alphabetical order). -```r +``` r nordic_cat <- factor( nordic_cat, levels = c( @@ -295,7 +295,7 @@ nordic_cat <- factor( nordic_cat ``` -```output +``` output [1] Norway Sweden Norway Denmark Sweden Levels: Norway Denmark Sweden ``` @@ -306,7 +306,7 @@ There is more than one way to reorder factors. Later in the lesson, we will use `fct_relevel()` function from `forcats` package to do the reordering. -```r +``` r library(forcats) nordic_cat <- fct_relevel( @@ -320,7 +320,7 @@ nordic_cat <- fct_relevel( nordic_cat ``` -```output +``` output [1] Norway Sweden Norway Denmark Sweden Levels: Norway Denmark Sweden ``` @@ -332,11 +332,11 @@ it shows the underlying values of each category. You can also see the structure in the environment tab of RStudio. -```r +``` r str(nordic_cat) ``` -```output +``` output Factor w/ 3 levels "Norway","Denmark",..: 1 3 1 2 3 ``` @@ -351,15 +351,15 @@ outside of this set, it will become an unknown/missing value detonated by -```r +``` r nordic_str ``` -```output +``` output [1] "Norway" "Sweden" "Norway" "Denmark" "Sweden" ``` -```r +``` r nordic_cat2 <- factor( nordic_str, levels = c("Norway", "Denmark") @@ -370,7 +370,7 @@ nordic_cat2 <- factor( nordic_cat2 ``` -```output +``` output [1] Norway Norway Denmark Levels: Norway Denmark ``` diff --git a/03-explore-data.md b/03-explore-data.md index bd4cdfa3..fa150a3a 100644 --- a/03-explore-data.md +++ b/03-explore-data.md @@ -54,7 +54,7 @@ For example, here is a figure depicting a data frame comprising a numeric, a cha We're gonna read in the `gapminder` data set with information about countries' size, GDP and average life expectancy in different years. -```r +``` r gapminder <- read.csv("data/gapminder_data.csv") ``` @@ -64,11 +64,11 @@ Let’s investigate the `gapminder` data frame a bit; the first thing we should It is important to see if all the variables (columns) have the data type that we require. For instance, a column might have numbers stored as characters, which would not allow us to make calculations with those numbers. -```r +``` r str(gapminder) ``` -```output +``` output 'data.frame': 1704 obs. of 6 variables: $ country : chr "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ... $ year : int 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ... @@ -88,11 +88,11 @@ There are multiple ways to explore a data set. Here are just a few examples: -```r +``` r head(gapminder) # shows first 6 rows of the data set ``` -```output +``` output country year pop continent lifeExp gdpPercap 1 Afghanistan 1952 8425333 Asia 28.801 779.4453 2 Afghanistan 1957 9240934 Asia 30.332 820.8530 @@ -102,11 +102,11 @@ head(gapminder) # shows first 6 rows of the data set 6 Afghanistan 1977 14880372 Asia 38.438 786.1134 ``` -```r +``` r summary(gapminder) # basic statistical information about each column. ``` -```output +``` output country year pop continent Length:1704 Min. :1952 Min. :6.001e+04 Length:1704 Class :character 1st Qu.:1966 1st Qu.:2.794e+06 Class :character @@ -123,21 +123,21 @@ summary(gapminder) # basic statistical information about each column. Max. :82.60 Max. :113523.1 ``` -```r +``` r # Information format differes by data type. nrow(gapminder) # returns number of rows in a dataset ``` -```output +``` output [1] 1704 ``` -```r +``` r ncol(gapminder) # returns number of columns in a dataset ``` -```output +``` output [1] 6 ``` @@ -147,7 +147,7 @@ When you're analyzing a data set, you often need to access its specific columns. One handy way to access a column is using it's name and a dollar sign `$`: -```r +``` r # This notation means: From dataset gapminder, give me column country. You can # see that the column accessed in this way is just a vector of characters. country_vec <- gapminder$country @@ -155,7 +155,7 @@ country_vec <- gapminder$country head(country_vec) ``` -```output +``` output [1] "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" [6] "Afghanistan" ``` @@ -170,13 +170,13 @@ Let's start manipulating the data. First, we will adapt our data set, by keeping only the columns we're interested in, using the `select()` function from the `dplyr` package: -```r +``` r year_country_gdp <- select(gapminder, year, country, gdpPercap) head(year_country_gdp) ``` -```output +``` output year country gdpPercap 1 1952 Afghanistan 779.4453 2 1957 Afghanistan 820.8530 @@ -196,14 +196,14 @@ In newer installation of `R` you can also find a notation `|>` . This pipe works The `select()` statement with pipe would look like that: -```r +``` r year_country_gdp <- gapminder %>% select(year, country, gdpPercap) head(year_country_gdp) ``` -```output +``` output year country gdpPercap 1 1952 Afghanistan 779.4453 2 1957 Afghanistan 820.8530 @@ -221,7 +221,7 @@ We already know how to select only the needed columns. But now, we also want to In the `gapminder` data set, we want to see the results from outside of Europe for the 21st century. -```r +``` r year_country_gdp_euro <- gapminder %>% filter(continent != "Europe" & year >= 2000) %>% select(year, country, gdpPercap) @@ -230,7 +230,7 @@ year_country_gdp_euro <- gapminder %>% head(year_country_gdp_euro) ``` -```output +``` output year country gdpPercap 1 2002 Afghanistan 726.7341 2 2007 Afghanistan 974.5803 @@ -258,7 +258,7 @@ year_country_gdp_eurasia <- gapminder %>% nrow(year_country_gdp_eurasia) ``` -```output +``` output [1] 756 ``` @@ -270,13 +270,13 @@ nrow(year_country_gdp_eurasia) So far, we have provided summary statistics on the whole dataset, selected columns, and filtered the observations. But often instead of doing that, we would like to know statistics about all of the continents, presented by group. -```r +``` r gapminder %>% # select the dataset group_by(continent) %>% # group by continent summarize(avg_gdpPercap = mean(gdpPercap)) # create basic stats ``` -```output +``` output # A tibble: 5 × 2 continent avg_gdpPercap @@ -306,7 +306,7 @@ gapminder %>% avg_lifeExp == max(avg_lifeExp)) ``` -```output +``` output # A tibble: 2 × 2 country avg_lifeExp @@ -322,13 +322,13 @@ gapminder %>% You can also group by multiple columns: -```r +``` r gapminder %>% group_by(continent, year) %>% summarize(avg_gdpPercap = mean(gdpPercap)) ``` -```output +``` output # A tibble: 60 × 3 # Groups: continent [5] continent year avg_gdpPercap @@ -348,7 +348,7 @@ gapminder %>% On top of this, you can also make multiple summaries of those groups: -```r +``` r gdp_pop_bycontinents_byyear <- gapminder %>% group_by(continent, year) %>% summarize( @@ -364,13 +364,13 @@ gdp_pop_bycontinents_byyear <- gapminder %>% If you need only a number of observations per group, you can use the `count()` function -```r +``` r gapminder %>% group_by(continent) %>% count() ``` -```output +``` output # A tibble: 5 × 2 # Groups: continent [5] continent n @@ -388,14 +388,14 @@ gapminder %>% Frequently you’ll want to create new columns based on the values in existing columns. For example, instead of only having the GDP per capita, we might want to create a new GDP variable and convert its units into Billions. For this, we’ll use `mutate()`. -```r +``` r gapminder_gdp <- gapminder %>% mutate(gdpBillion = gdpPercap * pop / 10^9) head(gapminder_gdp) ``` -```output +``` output country year pop continent lifeExp gdpPercap gdpBillion 1 Afghanistan 1952 8425333 Asia 28.801 779.4453 6.567086 2 Afghanistan 1957 9240934 Asia 30.332 820.8530 7.585449 diff --git a/04-intro-to-visualisation.md b/04-intro-to-visualisation.md index d8806a75..5e211d51 100644 --- a/04-intro-to-visualisation.md +++ b/04-intro-to-visualisation.md @@ -43,7 +43,7 @@ A fun part about `ggplot2` is that you can add layers to the plot to provide mor In the following parts of this workshop, you will use this package to visualize geospatial data. First, make sure that you have the following packages loaded. -```r +``` r library(tidyverse) library(terra) ``` @@ -51,7 +51,7 @@ library(terra) Now, lets plot the distribution of life expectancy in the `gapminder` dataset: -```r +``` r ggplot( data = gapminder, # data aes(x = lifeExp) # aesthetics layer @@ -69,7 +69,7 @@ already learned: `%>%` ( or `|>`) and follow them by ggplot syntax. Let's create another plot, this time only on a subset of observations: -```r +``` r gapminder %>% # we select a data set filter(year == 2007 & continent == "Americas") %>% # filter to keep one year and one continent ggplot(aes(x = country, y = gdpPercap)) + # the x and y axes represent values of columns @@ -82,7 +82,7 @@ Now, you can iteratively improve how the plot looks like. For example, you might want to flip it, to better display the labels. -```r +``` r gapminder %>% filter( year == 2007, @@ -104,7 +104,7 @@ Now the order of the levels will depend on another variable - GDP per capita. -```r +``` r gapminder %>% filter( year == 2007, @@ -122,7 +122,7 @@ Let's make things more colourful - let's represent the average life expectancy of a country by colour -```r +``` r gapminder %>% filter( year == 2007, @@ -145,7 +145,7 @@ readability and colorblind-proofness are the palettes available in the `viridis` package. -```r +``` r gapminder %>% filter( year == 2007, @@ -164,7 +164,7 @@ Maybe we don't need that much information about the life expectancy. We only want to know if it's below or above average. We will make use of the `if_else()` function inside `mutate()` to create a new column `lifeExpCat` with the value `high` if life expectancy is above average and `low` otherwise. Note the usage of the `if_else()` function: `if_else(, , )`. -```r +``` r p <- # this time let's save the plot in an object gapminder %>% filter(year == 2007 & @@ -193,7 +193,7 @@ like with any other object in `R`, if you want to see it, you need to call it. -```r +``` r p ``` @@ -204,7 +204,7 @@ Now we can make use of the saved object and add things to it. Let's also give it a title and name the axes: -```r +``` r p <- p + ggtitle("GDP per capita in Americas", subtitle = "Year 2007") + xlab("Country") + @@ -224,7 +224,7 @@ Once we are happy with our plot we can save it in a format of our choice. Remember to save it in the dedicated folder. -```r +``` r ggsave( plot = p, filename = here("fig_output", "plot_americas_2007.pdf") @@ -240,7 +240,7 @@ exploring the help documentation. We can do that by writing directly in the console: -```r +``` r ?ggsave ``` @@ -256,7 +256,7 @@ your analysis, you can then load directly that data set. Let's say we want to save the data only for Americas: -```r +``` r gapminder_amr_2007 <- gapminder %>% filter(year == 2007 & continent == "Americas") %>% mutate( diff --git a/09-open-and-plot-vector-layers.md b/09-open-and-plot-vector-layers.md index 9c66fa32..5bd9c477 100644 --- a/09-open-and-plot-vector-layers.md +++ b/09-open-and-plot-vector-layers.md @@ -41,7 +41,7 @@ In this lesson you will work with the `sf` package. Note that the `sf` package h First we need to load the packages we will use in this lesson. We will use the `tidyverse` package with which you are already familiar from the previous lesson. In addition, we need to load the [`sf`](https://r-spatial.github.io/sf/) package for working with spatial vector data. -```r +``` r library(tidyverse) # wrangle, reshape and visualize data library(sf) # work with spatial vector data ``` @@ -59,7 +59,7 @@ library(sf) # work with spatial vector data Let's start by opening a shapefile. Shapefiles are a common file format to store spatial vector data used in GIS software. Note that a shapefile consists of multiple files and it is important to keep them all together in the same location. We will read a shapefile with the administrative boundary of Delft with the function `st_read()` from the `sf` package. -```r +``` r boundary_Delft <- st_read("data/delft-boundary.shp", quiet = TRUE) ``` @@ -81,11 +81,11 @@ By default (with `quiet = FALSE`), the `st_read()` function provides a message w The `st_geometry_type()` function, for instance, gives us information about the geometry type, which in this case is `POLYGON`. -```r +``` r st_geometry_type(boundary_Delft) ``` -```output +``` output [1] POLYGON 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE ``` @@ -93,11 +93,11 @@ st_geometry_type(boundary_Delft) The `st_crs()` function returns the coordinate reference system (CRS) used by the shapefile, which in this case is `WGS 84` and has the unique reference code `EPSG: 4326`. -```r +``` r st_crs(boundary_Delft) ``` -```output +``` output Coordinate Reference System: User input: WGS 84 wkt: @@ -129,11 +129,11 @@ The `st_bbox()` function shows the extent of the layer. As `WGS 84` is a **geogr -```r +``` r st_bbox(boundary_Delft) ``` -```output +``` output xmin ymin xmax ymax 4.320218 51.966316 4.407911 52.032599 ``` @@ -142,51 +142,51 @@ We need a **projected CRS**, which in the case of the Netherlands is typically t To check the EPSG code of any CRS, we can check this website: https://epsg.io/ -```r +``` r boundary_Delft <- st_transform(boundary_Delft, 28992) st_crs(boundary_Delft)$Name ``` -```output +``` output [1] "Amersfoort / RD New" ``` -```r +``` r st_crs(boundary_Delft)$epsg ``` -```output +``` output [1] 28992 ``` Notice that the bounding box is measured in meters after the transformation. -```r +``` r st_bbox(boundary_Delft) ``` -```output +``` output xmin ymin xmax ymax 81743.00 442446.21 87703.78 449847.95 ``` -```r +``` r st_crs(boundary_Delft)$units_gdal ``` -```output +``` output [1] "metre" ``` We confirm the transformation by examining the reprojected shapefile. -```r +``` r boundary_Delft ``` -```output +``` output Simple feature collection with 1 feature and 1 field Geometry type: POLYGON Dimension: XY @@ -209,7 +209,7 @@ Read more about Coordinate Reference Systems [here](https://datacarpentry.org/or Now, let's plot this shapefile. You are already familiar with the `ggplot2` package from [Introduction to Visualisation](../episodes/04-intro-to-visualisation.Rmd). `ggplot2` has special `geom_` functions for spatial data. We will use the `geom_sf()` function for `sf` data. -```r +``` r ggplot(data = boundary_Delft) + geom_sf(size = 3, color = "black", fill = "cyan1") + labs(title = "Delft Administrative Boundary") + @@ -232,7 +232,7 @@ Read in `delft-streets.shp` and `delft-leisure.shp` and assign them to `lines_De :::::::::::::::::::::::: solution -```r +``` r lines_Delft <- st_read("data/delft-streets.shp") points_Delft <- st_read("data/delft-leisure.shp") ``` @@ -240,20 +240,20 @@ points_Delft <- st_read("data/delft-leisure.shp") We can check the type of type of geometry with the `st_geometry_type()` function. `lines_Delft` contains `"LINESTRING"` geometry and `points_Delft` is made of `"POINT"` geometries. -```r +``` r st_geometry_type(lines_Delft)[1] ``` -```output +``` output [1] LINESTRING 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE ``` -```r +``` r st_geometry_type(points_Delft)[2] ``` -```output +``` output [1] POINT 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE ``` @@ -261,39 +261,39 @@ st_geometry_type(points_Delft)[2] `lines_Delft` and `points_Delft` are in the correct CRS. -```r +``` r st_crs(lines_Delft)$epsg ``` -```output +``` output [1] 28992 ``` -```r +``` r st_crs(points_Delft)$epsg ``` -```output +``` output [1] 28992 ``` When looking at the bounding boxes with the `st_bbox()` function, we see the spatial extent of the two objects in a projected CRS using meters as units. `lines_Delft()` and `points_Delft` have similar extents. -```r +``` r st_bbox(lines_Delft) ``` -```output +``` output xmin ymin xmax ymax 81759.58 441223.13 89081.41 449845.81 ``` -```r +``` r st_bbox(points_Delft) ``` -```output +``` output xmin ymin xmax ymax 81863.21 442621.15 87370.15 449345.08 ``` diff --git a/10-explore-and-plot-by-vector-layer-attributes.md b/10-explore-and-plot-by-vector-layer-attributes.md index 36fdb584..6e86fd8e 100644 --- a/10-explore-and-plot-by-vector-layer-attributes.md +++ b/10-explore-and-plot-by-vector-layer-attributes.md @@ -30,11 +30,11 @@ After completing this episode, participants should be able to… Let's have a look at the content of the loaded data, starting with `lines_Delft`. In essence, an `"sf"` object is a data.frame with a "sticky" geometry column and some extra metadata, like the CRS, extent and geometry type we examined earlier. -```r +``` r lines_Delft ``` -```output +``` output Simple feature collection with 11244 features and 2 fields Geometry type: LINESTRING Dimension: XY @@ -57,22 +57,22 @@ First 10 features: This means that we can examine and manipulate them as data frames. For instance, we can look at the number of variables (columns in a data frame) with `ncol()`. -```r +``` r ncol(lines_Delft) ``` -```output +``` output [1] 3 ``` In the case of `point_Delft` those columns are `"osm_id"`, `"highway"` and `"geometry"`. We can check the names of the columns with the base R function `names()`. -```r +``` r names(lines_Delft) ``` -```output +``` output [1] "osm_id" "highway" "geometry" ``` @@ -88,11 +88,11 @@ attribute table. We can also preview the content of the object by looking at the first 6 rows with the `head()` function, which in the case of an `sf` object is similar to examining the object directly. -```r +``` r head (lines_Delft) ``` -```output +``` output Simple feature collection with 6 features and 2 fields Geometry type: LINESTRING Dimension: XY @@ -115,11 +115,11 @@ Using the `$` operator, we can examine the content of a single field of our line We can see the contents of the `highway` field of our lines feature: -```r +``` r head(lines_Delft$highway, 10) ``` -```output +``` output [1] "cycleway" "cycleway" "cycleway" "footway" "footway" "footway" [7] "service" "steps" "footway" "footway" ``` @@ -127,11 +127,11 @@ head(lines_Delft$highway, 10) To see only unique values within the `highway` field, we can use the `unique()` function. This function extracts all possible values of a character variable. -```r +``` r unique(lines_Delft$highway) ``` -```output +``` output [1] "cycleway" "footway" "service" "steps" [5] "residential" "unclassified" "construction" "secondary" [9] "busway" "living_street" "motorway_link" "tertiary" @@ -146,11 +146,11 @@ unique(lines_Delft$highway) R is also able to handle categorical variables called factors. With factors, we can use the `levels()` function to show unique values. To examine unique values of the `highway` variable this way, we have to first transform it into a factor with the `factor()` function: -```r +``` r levels(factor(lines_Delft$highway)) ``` -```output +``` output [1] "bridleway" "busway" "construction" "cycleway" [5] "footway" "living_street" "motorway" "motorway_link" [9] "path" "pedestrian" "platform" "primary" @@ -181,11 +181,11 @@ Explore the attributes associated with the `point_Delft` spatial object. 1. To find the number of attributes, we use the `ncol()` function: -```r +``` r ncol(point_Delft) ``` -```output +``` output [1] 3 ``` @@ -194,11 +194,11 @@ ncol(point_Delft) Using the `head()` function which displays 6 rows by default, we only see two values and `NA`s. -```r +``` r head(point_Delft) ``` -```output +``` output Simple feature collection with 6 features and 2 fields Geometry type: POINT Dimension: XY @@ -216,11 +216,11 @@ Projected CRS: Amersfoort / RD New We can increase the number of rows with the n argument (e.g., `head(n = 10)` to show 10 rows) until we see at least three distinct values in the leisure column. Note that printing an `sf` object will also display the first 10 rows. -```r +``` r head(point_Delft, 10) ``` -```output +``` output Simple feature collection with 10 features and 2 fields Geometry type: POINT Dimension: XY @@ -239,18 +239,18 @@ Projected CRS: Amersfoort / RD New 10 1148515039 playground POINT (84661.3 446818) ``` -```r +``` r # you might be lucky to see three distinct values ``` We have our answer (`sports_centre` is the third value), but in general this is not a good approach as the first rows might still have many `NA`s and three distinct values might still not be present in the first `n` rows of the data frame. To remove `NA`s, we can use the function `na.omit()` on the leisure column to remove `NA`s completely. Note that we use the `$` operator to examine the content of a single variable. -```r +``` r head(na.omit(point_Delft$leisure)) # this is better ``` -```output +``` output [1] "picnic_table" "marina" "marina" "sports_centre" [5] "sports_centre" "playground" ``` @@ -258,26 +258,26 @@ head(na.omit(point_Delft$leisure)) # this is better To show only unique values, we can use the `levels()` function on a factor to only see the first occurrence of each distinct value. Note `NA`s are dropped in this case and that we get the first three of the unique alphabetically ordered values. -```r +``` r head(levels(factor(point_Delft$leisure)), n = 3) ``` -```output +``` output [1] "dance" "dog_park" "escape_game" ``` -```r +``` r # this is even better ``` 3. To see a list of all attribute names, we can use the `names()` function. -```r +``` r names(point_Delft) ``` -```output +``` output [1] "osm_id" "leisure" "geometry" ``` @@ -291,7 +291,7 @@ names(point_Delft) We can use the `filter()` function to select a subset of features from a spatial object, just like with data frames. Let's select only cycleways from our street data. -```r +``` r cycleway_Delft <- lines_Delft %>% filter(highway == "cycleway") ``` @@ -299,26 +299,26 @@ cycleway_Delft <- lines_Delft %>% Our subsetting operation reduces the number of features from 11244 to 1397. -```r +``` r nrow(lines_Delft) ``` -```output +``` output [1] 11244 ``` -```r +``` r nrow(cycleway_Delft) ``` -```output +``` output [1] 1397 ``` This can be useful, for instance, to calculate the total length of cycleways. -```r +``` r cycleway_Delft <- cycleway_Delft %>% mutate(length = st_length(.)) @@ -326,7 +326,7 @@ cycleway_Delft %>% summarise(total_length = sum(length)) ``` -```output +``` output Simple feature collection with 1 feature and 1 field Geometry type: MULTILINESTRING Dimension: XY @@ -339,7 +339,7 @@ Projected CRS: Amersfoort / RD New Now we can plot only the cycleways. -```r +``` r ggplot(data = cycleway_Delft) + geom_sf() + labs( @@ -370,11 +370,11 @@ Challenge: Now with motorways (and pedestrian streets) 1. To create the new object, we first need to see which value of the `highway` column holds motorways. There is a value called `motorway`. -```r +``` r unique(lines_Delft$highway) ``` -```output +``` output [1] "cycleway" "footway" "service" "steps" [5] "residential" "unclassified" "construction" "secondary" [9] "busway" "living_street" "motorway_link" "tertiary" @@ -386,14 +386,14 @@ unique(lines_Delft$highway) We extract only the features with the value `motorway`. -```r +``` r motorway_Delft <- lines_Delft %>% filter(highway == "motorway") motorway_Delft ``` -```output +``` output Simple feature collection with 48 features and 2 fields Geometry type: LINESTRING Dimension: XY @@ -416,7 +416,7 @@ First 10 features: 2. There are 48 features with the value `motorway`. -```r +``` r motorway_Delft_length <- motorway_Delft %>% mutate(length = st_length(.)) %>% select(everything(), geometry) %>% @@ -426,18 +426,18 @@ motorway_Delft_length <- motorway_Delft %>% 3. The total length of motorways is 14877.4361477941. -```r +``` r nrow(motorway_Delft) ``` -```output +``` output [1] 48 ``` 4. Plot the motorways. -```r +``` r ggplot(data = motorway_Delft) + geom_sf(linewidth = 1.5) + labs( @@ -452,7 +452,7 @@ ggplot(data = motorway_Delft) + 5. Follow the same steps with pedestrian streets. -```r +``` r pedestrian_Delft <- lines_Delft %>% filter(highway == "pedestrian") @@ -462,7 +462,7 @@ pedestrian_Delft %>% summarise(total_length = sum(length)) ``` -```output +``` output Simple feature collection with 1 feature and 1 field Geometry type: MULTILINESTRING Dimension: XY @@ -472,16 +472,16 @@ Projected CRS: Amersfoort / RD New 1 12778.84 [m] MULTILINESTRING ((85538.03 ... ``` -```r +``` r nrow(pedestrian_Delft) ``` -```output +``` output [1] 234 ``` -```r +``` r ggplot() + geom_sf(data = pedestrian_Delft) + labs( @@ -502,11 +502,11 @@ ggplot() + Let's say that we want to color different road types with different colors and that we want to determine those colors. -```r +``` r unique(lines_Delft$highway) ``` -```output +``` output [1] "cycleway" "footway" "service" "steps" [5] "residential" "unclassified" "construction" "secondary" [9] "busway" "living_street" "motorway_link" "tertiary" @@ -519,7 +519,7 @@ unique(lines_Delft$highway) If we look at all the unique values of the highway field of our street network we see more than 20 values. Let's focus on a subset of four values to illustrate the use of distinct colors. We use a piped expression in which we only filter the rows of our data frame that have one of the four given values `"motorway"`, `"primary"`, `"secondary"`, and `"cycleway"`. Note that we do this with the `%in` operator which is a more compact equivalent of a series of conditions joined by the `|` (or) operator. We also make sure that the highway column is a factor column. -```r +``` r road_types <- c("motorway", "primary", "secondary", "cycleway") lines_Delft_selection <- lines_Delft %>% @@ -530,14 +530,14 @@ lines_Delft_selection <- lines_Delft %>% Next we define the four colors we want to use, one for each type of road in our vector object. Note that in R you can use named colors like `"blue"`, `"green"`, `"navy"`, and `"purple"`. If you are using RStudio, you will see the named colors previewed in line. A full list of named colors can be listed with the `colors()` function. -```r +``` r road_colors <- c("blue", "green", "navy", "purple") ``` We can use the defined color palette in ggplot. -```r +``` r ggplot(data = lines_Delft_selection) + geom_sf(aes(color = highway)) + scale_color_manual(values = road_colors) + @@ -556,12 +556,12 @@ ggplot(data = lines_Delft_selection) + Earlier we adjusted the line width universally. We can also adjust line widths for every factor level. Note that in this case the `size` argument, like the `color` argument, are within the `aes()` mapping function. This means that the values of that visual property will be mapped from a variable of the object that is being plotted. -```r +``` r line_widths <- c(1, 0.75, 0.5, 0.25) ``` -```r +``` r ggplot(data = lines_Delft_selection) + geom_sf(aes(color = highway, linewidth = highway)) + scale_color_manual(values = road_colors) + @@ -593,21 +593,21 @@ Let’s create another plot where we show the different line types with the foll ::: solution -```r +``` r levels(factor(lines_Delft_selection$highway)) ``` -```output +``` output [1] "motorway" "primary" "secondary" "cycleway" ``` -```r +``` r line_width <- c(0.25, 0.75, 0.5, 1) ``` -```r +``` r ggplot(data = lines_Delft_selection) + geom_sf(aes(linewidth = highway)) + scale_linewidth_manual(values = line_width) + @@ -630,7 +630,7 @@ ggplot(data = lines_Delft_selection) + Let’s add a legend to our plot. We will use the `road_colors` object that we created above to color the legend. We can customize the appearance of our legend by manually setting different parameters. -```r +``` r p1 <- ggplot(data = lines_Delft_selection) + geom_sf(aes(color = highway), linewidth = 1.5) + scale_color_manual(values = road_colors) + @@ -651,7 +651,7 @@ p1 -```r +``` r p2 <- p1 + theme( legend.text = element_text(size = 20), @@ -677,20 +677,20 @@ Create a plot that emphasizes only roads where bicycles are allowed. To emphasiz ::: solution -```r +``` r class(lines_Delft_selection$highway) ``` -```output +``` output [1] "factor" ``` -```r +``` r levels(factor(lines_Delft$highway)) ``` -```output +``` output [1] "bridleway" "busway" "construction" "cycleway" [5] "footway" "living_street" "motorway" "motorway_link" [9] "path" "pedestrian" "platform" "primary" @@ -701,7 +701,7 @@ levels(factor(lines_Delft$highway)) ``` -```r +``` r # First, create a data frame with only roads where bicycles # are allowed lines_Delft_bicycle <- lines_Delft %>% @@ -739,11 +739,11 @@ Create a map of the municipal boundaries in the Netherlands using the data locat ::: solution -```r +``` r municipal_boundaries_NL <- st_read("data/nl-gemeenten.shp") ``` -```output +``` output Reading layer `nl-gemeenten' from data source `/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/nl-gemeenten.shp' using driver `ESRI Shapefile' @@ -755,11 +755,11 @@ Projected CRS: Amersfoort / RD New ``` -```r +``` r str(municipal_boundaries_NL) ``` -```output +``` output Classes 'sf' and 'data.frame': 344 obs. of 7 variables: $ identifica: chr "GM0014" "GM0034" "GM0037" "GM0047" ... $ naam : chr "Groningen" "Almere" "Stadskanaal" "Veendam" ... @@ -776,18 +776,18 @@ Classes 'sf' and 'data.frame': 344 obs. of 7 variables: ..- attr(*, "names")= chr [1:6] "identifica" "naam" "code" "ligtInProv" ... ``` -```r +``` r levels(factor(municipal_boundaries_NL$ligtInPr_1)) ``` -```output +``` output [1] "Drenthe" "Flevoland" "Fryslân" "Gelderland" [5] "Groningen" "Limburg" "Noord-Brabant" "Noord-Holland" [9] "Overijssel" "Utrecht" "Zeeland" "Zuid-Holland" ``` -```r +``` r ggplot(data = municipal_boundaries_NL) + geom_sf(aes(color = ligtInPr_1), linewidth = 1) + labs(title = "Contiguous NL Municipal Boundaries") + diff --git a/11-plot-multiple-shape-files.md b/11-plot-multiple-shape-files.md index 314f01f7..44315d91 100644 --- a/11-plot-multiple-shape-files.md +++ b/11-plot-multiple-shape-files.md @@ -41,7 +41,7 @@ We will create a plot that combines our leisure locations (`point_Delft`), munic To begin, we will create a plot with the site boundary as the first layer. Then layer the leisure locations and street data on top using `+`. -```r +``` r ggplot() + geom_sf( data = boundary_Delft, @@ -63,7 +63,7 @@ ggplot() + Next, let’s build a custom legend using the functions `scale_color_manual()` and `scale_fill_manual()`. -```r +``` r leisure_colors <- rainbow(15) point_Delft$leisure <- factor(point_Delft$leisure) @@ -98,7 +98,7 @@ ggplot() + -```r +``` r ggplot() + geom_sf( data = boundary_Delft, @@ -143,12 +143,12 @@ Modify the plot above. Tell R to plot each point, using a different symbol of sh :::::::::::::::::::::::: solution -```r +``` r leisure_locations_selection <- st_read("data/delft-leisure.shp") %>% filter(leisure %in% c("playground", "picnic_table")) ``` -```output +``` output Reading layer `delft-leisure' from data source `/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/delft-leisure.shp' using driver `ESRI Shapefile' @@ -160,21 +160,21 @@ Projected CRS: Amersfoort / RD New ``` -```r +``` r levels(factor(leisure_locations_selection$leisure)) ``` -```output +``` output [1] "picnic_table" "playground" ``` -```r +``` r blue_orange <- c("cornflowerblue", "darkorange") ``` -```r +``` r ggplot() + geom_sf( data = lines_Delft_selection, @@ -209,7 +209,7 @@ ggplot() + -```r +``` r ggplot() + geom_sf( data = lines_Delft_selection, diff --git a/12-handling-spatial-projection-and-crs.md b/12-handling-spatial-projection-and-crs.md index d2c6ec71..42b8cc07 100644 --- a/12-handling-spatial-projection-and-crs.md +++ b/12-handling-spatial-projection-and-crs.md @@ -24,11 +24,11 @@ After completing this episode, participants should be able to… ## Working with spatial data from different sources -```r +``` r municipal_boundary_NL <- st_read("data/nl-gemeenten.shp") ``` -```output +``` output Reading layer `nl-gemeenten' from data source `/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/nl-gemeenten.shp' using driver `ESRI Shapefile' @@ -40,7 +40,7 @@ Projected CRS: Amersfoort / RD New ``` -```r +``` r ggplot() + geom_sf(data = municipal_boundary_NL) + labs(title = "Map of Contiguous NL Municipal Boundaries") + @@ -52,11 +52,11 @@ ggplot() + We can add a country boundary layer to make it look nicer. If we specify a thicker line width using size = 2 for the country boundary layer, it will make our map pop! -```r +``` r country_boundary_NL <- st_read("data/nl-boundary.shp") ``` -```output +``` output Reading layer `nl-boundary' from data source `/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/nl-boundary.shp' using driver `ESRI Shapefile' @@ -68,7 +68,7 @@ Projected CRS: Amersfoort / RD New ``` -```r +``` r ggplot() + geom_sf( data = country_boundary_NL, @@ -86,34 +86,34 @@ ggplot() + -```r +``` r # st_crs(point_Delft) ``` -```r +``` r st_crs(municipal_boundary_NL)$epsg ``` -```output +``` output [1] 28992 ``` -```r +``` r st_crs(country_boundary_NL)$epsg ``` -```output +``` output [1] 28992 ``` -```r +``` r boundary_Delft <- st_read("data/delft-boundary.shp") ``` -```output +``` output Reading layer `delft-boundary' from data source `/home/runner/work/r-geospatial-urban/r-geospatial-urban/site/built/data/delft-boundary.shp' using driver `ESRI Shapefile' @@ -124,20 +124,20 @@ Bounding box: xmin: 4.320218 ymin: 51.96632 xmax: 4.407911 ymax: 52.0326 Geodetic CRS: WGS 84 ``` -```r +``` r st_crs(boundary_Delft)$epsg ``` -```output +``` output [1] 4326 ``` -```r +``` r boundary_Delft <- st_transform(boundary_Delft, 28992) ``` -```r +``` r ggplot() + geom_sf( data = country_boundary_NL, @@ -175,13 +175,13 @@ Create a map of South Holland as follows: ::: solution -```r +``` r boundary_ZH <- municipal_boundary_NL %>% filter(ligtInPr_1 == "Zuid-Holland") ``` -```r +``` r ggplot() + geom_sf( data = boundary_ZH, @@ -222,7 +222,7 @@ ggplot() + To save a file, use the `st_write()` function from the `sf` package. Although `sf` guesses the driver needed for a specified output file name from its extension, this can be made explicitly via the `driver` argument. In our case `driver = "ESRI Shapefile"` ensures that the output is correctly saved as a `.shp` file. -```r +``` r st_write(leisure_locations_selection, "data/leisure_locations_selection.shp", driver = "ESRI Shapefile" diff --git a/13-intro-to-raster-data.md b/13-intro-to-raster-data.md index c66e492a..aec7f0ee 100644 --- a/13-intro-to-raster-data.md +++ b/13-intro-to-raster-data.md @@ -37,7 +37,7 @@ In this lesson, we will work with raster data. We will start with an introductio We continue to work with the `tidyverse` package and we will use the `terra` package to work with raster data. Make sure that you have those packages loaded. -```r +``` r library(tidyverse) library(terra) ``` @@ -58,11 +58,11 @@ In this and lesson, we will use: We will be working with a series of GeoTIFF files in this lesson. The GeoTIFF format contains a set of embedded tags with metadata about the raster data. We can use the function `describe()` from the `terra` package to get information about our raster data before we read that data into R. It is recommended to do this before importing your data. We first examine the file `tud-dsm-5m.tif`. -```r +``` r describe("data/tud-dsm-5m.tif") ``` -```output +``` output [1] "Driver: GTiff/GeoTIFF" [2] "Files: data/tud-dsm-5m.tif" [3] "Size is 722, 386" @@ -134,12 +134,12 @@ To improve code readability, use file and object names that make it clear what i First we will load our raster file into R and view the data structure. -```r +``` r DSM_TUD <- rast("data/tud-dsm-5m.tif") DSM_TUD ``` -```output +``` output class : SpatRaster dimensions : 386, 722, 1 (nrow, ncol, nlyr) resolution : 5, 5 (x, y) @@ -151,15 +151,15 @@ name : tud-dsm-5m The information above includes a report on dimension, resolution, extent and CRS, but no information about the values. Similar to other data structures in R like vectors and data frame columns, descriptive statistics for raster data can be retrieved with the `summary()` function. -```r +``` r summary(DSM_TUD) ``` -```warning +``` warning Warning: [summary] used a sample ``` -```output +``` output tud.dsm.5m Min. :-5.2235 1st Qu.:-0.7007 @@ -172,11 +172,11 @@ Warning: [summary] used a sample This output gives us information about the range of values in the DSM. We can see, for instance, that the lowest elevation is `-5.2235`, the highest is `89.7838`. But note the warning. Unless you force R to calculate these statistics using every cell in the raster, it will take a random sample of 100,000 cells and calculate from them instead. To force calculation all the values, you can use the function `values`: -```r +``` r summary(values(DSM_TUD)) ``` -```output +``` output tud-dsm-5m Min. :-5.3907 1st Qu.:-0.7008 @@ -191,18 +191,18 @@ With a summary on all cells of the raster, the values range from a smaller minim To visualise the DSM in R using `ggplot2`, we need to convert it to a data frame. We learned about data frames in an [earlier lesson](../episodes/03-explore-data.Rmd). The `terra` package has the built-in function `as.data.frame()` for conversion to a data frame. -```r +``` r DSM_TUD_df <- as.data.frame(DSM_TUD, xy = TRUE) ``` Now when we view the structure of our data, we will see a standard data frame format. -```r +``` r str(DSM_TUD_df) ``` -```output +``` output 'data.frame': 278692 obs. of 3 variables: $ x : num 83568 83572 83578 83582 83588 ... $ y : num 447178 447178 447178 447178 447178 ... @@ -212,7 +212,7 @@ str(DSM_TUD_df) We can use `ggplot()` to plot this data with a specific `geom_` function called `geom_raster()`. We will make the colour scale in our plot colour-blindness friendly with `scale_fill_viridis_c`, introduced in an [earlier lesson](../episodes/04-intro-to-visualisation.Rmd). We will also use the `coord_quickmap()` function to use an approximate Mercator projection for our plots. This approximation is suitable for small areas that are not too close to the poles. Other coordinate systems are available in `ggplot2` if needed, you can learn about them at their help page `?coord_map`. -```r +``` r ggplot() + geom_raster(data = DSM_TUD_df , aes(x = x, y = y, fill = `tud-dsm-5m`)) + scale_fill_viridis_c(option = "H") + # `option = "H"` provides a contrasting colour scale @@ -245,11 +245,11 @@ Now we will see how features of the CRS appear in our data file and what meaning We can view the CRS string associated with our R object using the `crs()` function. -```r +``` r crs(DSM_TUD, proj = TRUE) ``` -```output +``` output [1] "+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs" ``` @@ -277,11 +277,11 @@ It is useful to know the minimum and maximum values of a raster dataset. In this Raster statistics are often calculated and embedded in a GeoTIFF for us. We can view these values: -```r +``` r minmax(DSM_TUD) ``` -```output +``` output tud-dsm-5m min Inf max -Inf @@ -293,26 +293,26 @@ max -Inf If the `min` and `max` values are `Inf` and `-Inf` respectively, it means that they haven't been calculated. We can calculate them using the `setMinMax()` function. -```r +``` r DSM_TUD <- setMinMax(DSM_TUD) ``` ::: -```r +``` r min(values(DSM_TUD)) ``` -```output +``` output [1] -5.39069 ``` -```r +``` r max(values(DSM_TUD)) ``` -```output +``` output [1] 92.08102 ``` @@ -330,11 +330,11 @@ The Digital Surface Model object (`DSM_TUD`) that we’ve been working with is a A raster dataset can contain one or more bands. We can view the number of bands in a raster using the `nlyr()` function. -```r +``` r nlyr(DSM_TUD) ``` -```output +``` output [1] 1 ``` This dataset has only 1 band. However, raster data can also be multi-band, meaning that one raster file contains data for more than one variable or time period for each cell. We will discuss multi-band raster data in a [later episode](../episodes/17-work-with-multi-band-rasters.Rmd). @@ -348,12 +348,12 @@ A histogram can be used to inspect the distribution of raster values visually. I We can inspect the distribution of values contained in a raster using the `ggplot2` function `geom_histogram()`. Histograms are often useful in identifying outliers and bad data values in our raster data. -```r +``` r ggplot() + geom_histogram(data = DSM_TUD_df, aes(`tud-dsm-5m`)) ``` -```output +``` output `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ``` @@ -368,7 +368,7 @@ Notice that a message is displayed when R creates the histogram: This message is caused by a default setting in `geom_histogram()` enforcing that there are 30 bins for the data. We can define the number of bins we want in the histogram by giving another value to the `bins` argument. 60 bins, for instance, will give a more detailed histogram of the same distribution: -```r +``` r ggplot() + geom_histogram(data = DSM_TUD_df, aes(`tud-dsm-5m`), bins = 60) ``` @@ -393,11 +393,11 @@ Note that this file is a hillshade raster. We will learn about hillshades in the ::: solution -```r +``` r describe("data/tud-dsm-5m-hill.tif") ``` -```output +``` output [1] "Driver: GTiff/GeoTIFF" [2] "Files: data/tud-dsm-5m-hill.tif" [3] "Size is 722, 386" diff --git a/14-plot-raster-data.md b/14-plot-raster-data.md index 168f9421..fe1a2ca8 100644 --- a/14-plot-raster-data.md +++ b/14-plot-raster-data.md @@ -37,7 +37,7 @@ In this part, we will plot our raster object using `ggplot2` with customized col In the previous plot, our DSM was coloured with a continuous colour range. For clarity and visibility, we may prefer to view the data “symbolized” or coloured according to ranges of values. This is comparable to a “classified” map. For that, we need to tell `ggplot()` how many groups to break our data into and where those breaks should be. To make these decisions, it is useful to first explore the distribution of the data using a bar plot. To begin with, we will use `dplyr`’s `mutate()` function combined with `cut()` to split the data into 3 bins. -```r +``` r DSM_TUD_df <- DSM_TUD_df %>% mutate(fct_elevation = cut(`tud-dsm-5m`, breaks = 3)) @@ -49,22 +49,22 @@ ggplot() + To see the cut-off values for the groups, we can ask for the levels of `fct_elevation`: -```r +``` r levels(DSM_TUD_df$fct_elevation) ``` -```output +``` output [1] "(-5.49,27.1]" "(27.1,59.6]" "(59.6,92.2]" ``` And we can get the count of values (that is, number of pixels) in each group using `dplyr`’s `count()` function: -```r +``` r DSM_TUD_df %>% count(fct_elevation) ``` -```output +``` output fct_elevation n 1 (-5.49,27.1] 277100 2 (27.1,59.6] 1469 @@ -74,7 +74,7 @@ DSM_TUD_df %>% We might prefer to customize the cut-off values for these groups. Lets round the cut-off values so that we have groups for the ranges of -10 - 0m, 0 - 5m, and 5 - 100m. To implement this we will give `cut()` a numeric vector of break points instead of the number of breaks we want. -```r +``` r custom_bins <- c(-10, 0, 5, 100) DSM_TUD_df <- DSM_TUD_df %>% @@ -83,7 +83,7 @@ DSM_TUD_df <- DSM_TUD_df %>% levels(DSM_TUD_df$fct_elevation_cb) ``` -```output +``` output [1] "(-10,0]" "(0,5]" "(5,100]" ``` @@ -99,7 +99,7 @@ The bin intervals are shown using `(` to mean exclusive and `]` to mean inclusiv And now we can plot our bar plot again, using the new groups: -```r +``` r ggplot() + geom_bar(data = DSM_TUD_df, aes(fct_elevation_cb)) ``` @@ -108,12 +108,12 @@ ggplot() + And we can get the count of values in each group in the same way we did before: -```r +``` r DSM_TUD_df %>% count(fct_elevation_cb) ``` -```output +``` output fct_elevation_cb n 1 (-10,0] 113877 2 (0,5] 101446 @@ -122,7 +122,7 @@ DSM_TUD_df %>% We can use those groups to plot our raster data, with each group being a different colour: -```r +``` r ggplot() + geom_raster(data = DSM_TUD_df , aes(x = x, y = y, fill = fct_elevation_cb)) + coord_quickmap() @@ -132,17 +132,17 @@ ggplot() + The plot above uses the default colours inside `ggplot2` for raster objects. We can specify our own colours to make the plot look a little nicer. R has a built in set of colours for plotting terrain, which are built in to the `terrain.colors()` function. Since we have three bins, we want to create a 3-colour palette: -```r +``` r terrain.colors(3) ``` -```output +``` output [1] "#00A600" "#ECB176" "#F2F2F2" ``` The `terrain.colors()` function returns hex colours - each of these character strings represents a colour. To use these in our map, we pass them across using the `scale_fill_manual()` function. -```r +``` r ggplot() + geom_raster(data = DSM_TUD_df , aes(x = x, y = y, fill = fct_elevation_cb)) + scale_fill_manual(values = terrain.colors(3)) + @@ -157,7 +157,7 @@ If we need to create multiple plots using the same colour palette, we can create We can give the legend a more meaningful title by passing a value to the `name` argument of the `scale_fill_manual()` function. -```r +``` r my_col <- terrain.colors(3) ggplot() + @@ -170,7 +170,7 @@ ggplot() + The axis labels x and y are not necessary, so we can turn them off by passing `element_blank()` to the relevant part of the `theme()` function. -```r +``` r ggplot() + geom_raster(data = DSM_TUD_df , aes(x = x, y = y, fill = fct_elevation_cb)) + @@ -194,19 +194,19 @@ Create a plot of the TU Delft Digital Surface Model (`DSM_TUD`) that has: ::: solution -```r +``` r DSM_TUD_df <- DSM_TUD_df %>% mutate(fct_elevation_6 = cut(`tud-dsm-5m`, breaks = 6)) levels(DSM_TUD_df$fct_elevation_6) ``` -```output +``` output [1] "(-5.49,10.9]" "(10.9,27.1]" "(27.1,43.3]" "(43.3,59.6]" "(59.6,75.8]" [6] "(75.8,92.2]" ``` -```r +``` r my_col <- terrain.colors(6) ggplot() + diff --git a/15-reproject-raster-data.md b/15-reproject-raster-data.md index 51140727..0a893a32 100644 --- a/15-reproject-raster-data.md +++ b/15-reproject-raster-data.md @@ -53,21 +53,21 @@ Read more about layering raster [here](https://datacarpentry.org/r-raster-vector First, we need to import the DTM and DTM hillshade data. -```r +``` r DTM_TUD <- rast("data/tud-dtm-5m.tif") DTM_hill_TUD <- rast("data/tud-dtm-5m-hill-WGS84.tif") ``` Next, we will convert each of these datasets to a data frame for plotting with `ggplot`. -```r +``` r DTM_TUD_df <- as.data.frame(DTM_TUD, xy = TRUE) DTM_hill_TUD_df <- as.data.frame(DTM_hill_TUD, xy = TRUE) ``` Now we can create a map of the DTM layered over the hillshade. -```r +``` r ggplot() + geom_raster(data = DTM_TUD_df , aes(x = x, y = y, @@ -83,7 +83,7 @@ ggplot() + Our results are curious - neither the Digital Terrain Model (`DTM_TUD_df`) nor the DTM Hillshade (`DTM_hill_TUD_df`) plotted. Let’s try to plot the DTM on its own to make sure there are data there. -```r +``` r ggplot() + geom_raster(data = DTM_TUD_df, aes(x = x, y = y, @@ -99,7 +99,7 @@ Our DTM seems to contain data and plots just fine. Next we plot the DTM Hillshade on its own to see whether everything is OK. -```r +``` r ggplot() + geom_raster(data = DTM_hill_TUD_df, aes(x = x, y = y, @@ -120,11 +120,11 @@ View the CRS for each of these two datasets. What projection does each use? ::: solution -```r +``` r crs(DTM_TUD, parse = TRUE) ``` -```output +``` output [1] "PROJCRS[\"Amersfoort / RD New\"," [2] " BASEGEOGCRS[\"Amersfoort\"," [3] " DATUM[\"Amersfoort\"," @@ -166,11 +166,11 @@ crs(DTM_TUD, parse = TRUE) ``` -```r +``` r crs(DTM_hill_TUD, parse = TRUE) ``` -```output +``` output [1] "GEOGCRS[\"WGS 84\"," [2] " DATUM[\"World Geodetic System 1984\"," [3] " ELLIPSOID[\"WGS 84\",6378137,298.257223563," @@ -217,18 +217,18 @@ We want the CRS of our hillshade to match the `DTM_TUD` raster. We can thus assi First we will reproject our `DTM_hill_TUD` raster data to match the `DTM_TUD` raster CRS: -```r +``` r DTM_hill_EPSG28992_TUD <- project(DTM_hill_TUD, crs(DTM_TUD)) ``` Now we can compare the CRS of our original DTM hillshade and our new DTM hillshade, to see how they are different. -```r +``` r crs(DTM_hill_EPSG28992_TUD, parse = TRUE) ``` -```output +``` output [1] "PROJCRS[\"Amersfoort / RD New\"," [2] " BASEGEOGCRS[\"Amersfoort\"," [3] " DATUM[\"Amersfoort\"," @@ -269,11 +269,11 @@ crs(DTM_hill_EPSG28992_TUD, parse = TRUE) [38] " ID[\"EPSG\",28992]]" ``` -```r +``` r crs(DTM_hill_TUD, parse = TRUE) ``` -```output +``` output [1] "GEOGCRS[\"WGS 84\"," [2] " DATUM[\"World Geodetic System 1984\"," [3] " ELLIPSOID[\"WGS 84\",6378137,298.257223563," @@ -292,19 +292,19 @@ crs(DTM_hill_TUD, parse = TRUE) We can also compare the extent of the two objects. -```r +``` r ext(DTM_hill_EPSG28992_TUD) ``` -```output +``` output SpatExtent : 83537.3768729672, 87200.5199626113, 445202.584641046, 447230.395994242 (xmin, xmax, ymin, ymax) ``` -```r +``` r ext(DTM_hill_TUD) ``` -```output +``` output SpatExtent : 4.34674898644234, 4.39970436836596, 51.9910492930106, 52.0088368700157 (xmin, xmax, ymin, ymax) ``` @@ -313,26 +313,26 @@ SpatExtent : 4.34674898644234, 4.39970436836596, 51.9910492930106, 52.0088368700 Let’s next have a look at the resolution of our reprojected hillshade versus our original data. -```r +``` r res(DTM_hill_EPSG28992_TUD) ``` -```output +``` output [1] 5.03179 5.03179 ``` -```r +``` r res(DTM_TUD) ``` -```output +``` output [1] 5 5 ``` These two resolutions are different, but they’re representing the same data. We can tell R to force our newly reprojected raster to be the same as `DTM_TUD` by adding a line of code `res = res(DTM_TUD)` within the `project()` function. -```r +``` r DTM_hill_EPSG28992_TUD <- project(DTM_hill_TUD, crs(DTM_TUD), res = res(DTM_TUD)) @@ -340,31 +340,31 @@ DTM_hill_EPSG28992_TUD <- project(DTM_hill_TUD, Now both our resolutions and our CRSs match, so we can plot these two data sets together. Let’s double-check our resolution to be sure: -```r +``` r res(DTM_hill_EPSG28992_TUD) ``` -```output +``` output [1] 5 5 ``` -```r +``` r res(DTM_TUD) ``` -```output +``` output [1] 5 5 ``` For plotting with `ggplot()`, we will need to create a data frame from our newly reprojected raster. -```r +``` r DTM_hill_TUD_2_df <- as.data.frame(DTM_hill_EPSG28992_TUD, xy = TRUE) ``` We can now create a plot of this data. -```r +``` r ggplot() + geom_raster(data = DTM_TUD_df , aes(x = x, y = y, diff --git a/16-raster-calculations.md b/16-raster-calculations.md index b4efc434..cc23f237 100644 --- a/16-raster-calculations.md +++ b/16-raster-calculations.md @@ -45,11 +45,11 @@ For this episode, we will use the DTM and DSM data which we already have loaded We use the `describe()` function to view information about the DTM and DSM data files. -```r +``` r describe("data/tud-dtm-5m.tif") ``` -```output +``` output [1] "Driver: GTiff/GeoTIFF" [2] "Files: data/tud-dtm-5m.tif" [3] "Size is 722, 386" @@ -108,11 +108,11 @@ describe("data/tud-dtm-5m.tif") [56] "Band 1 Block=722x2 Type=Float32, ColorInterp=Gray" ``` -```r +``` r describe("data/tud-dsm-5m.tif") ``` -```output +``` output [1] "Driver: GTiff/GeoTIFF" [2] "Files: data/tud-dsm-5m.tif" [3] "Size is 722, 386" @@ -173,7 +173,7 @@ describe("data/tud-dsm-5m.tif") We’ve already loaded and worked with these two data files in earlier episodes. Let’s plot them each once more to remind ourselves what this data looks like. First we’ll plot the DTM elevation data: -```r +``` r ggplot() + geom_raster(data = DTM_TUD_df , aes(x = x, y = y, fill = `tud-dtm-5m`)) + @@ -185,7 +185,7 @@ We’ve already loaded and worked with these two data files in earlier episodes. And then the DSM elevation data: -```r +``` r ggplot() + geom_raster(data = DSM_TUD_df , aes(x = x, y = y, fill = `tud-dsm-5m`)) + @@ -203,14 +203,14 @@ Let’s subtract the DTM from the DSM to create a Canopy Height Model. After sub -```r +``` r CHM_TUD <- DSM_TUD - DTM_TUD CHM_TUD_df <- as.data.frame(CHM_TUD, xy = TRUE) ``` We can now plot the output CHM. -```r +``` r ggplot() + geom_raster(data = CHM_TUD_df , aes(x = x, y = y, fill = `tud-dsm-5m`)) + @@ -222,7 +222,7 @@ We can now plot the output CHM. Let’s have a look at the distribution of values in our newly created Canopy Height Model (CHM). -```r +``` r ggplot(CHM_TUD_df) + geom_histogram(aes(`tud-dsm-5m`)) ``` @@ -243,30 +243,30 @@ It’s often a good idea to explore the range of values in a raster dataset just ::: solution -```r +``` r min(CHM_TUD_df$`tud-dsm-5m`, na.rm = TRUE) ``` -```output +``` output [1] -3.638057 ``` -```r +``` r max(CHM_TUD_df$`tud-dsm-5m`, na.rm = TRUE) ``` -```output +``` output [1] 92.08102 ``` -```r +``` r ggplot(CHM_TUD_df) + geom_histogram(aes(`tud-dsm-5m`)) ``` -```r +``` r custom_bins <- c(-5, 0, 10, 20, 30, 100) CHM_TUD_df <- CHM_TUD_df %>% mutate(canopy_discrete = cut(`tud-dsm-5m`, breaks = custom_bins)) @@ -308,7 +308,7 @@ When we write this raster object to a GeoTIFF file we’ll name it `CHM_TUD.tiff We will specify the output format (“GTiff”), the no data value `NAflag = -9999`. We will also tell R to overwrite any data that is already in a file of the same name. -```r +``` r writeRaster(CHM_TUD, "fig/CHM_TUD.tiff", filetype="GTiff", overwrite=TRUE, diff --git a/17-work-with-multi-band-rasters.md b/17-work-with-multi-band-rasters.md index 854d3bf4..a275ee88 100644 --- a/17-work-with-multi-band-rasters.md +++ b/17-work-with-multi-band-rasters.md @@ -42,13 +42,13 @@ In this episode, the multi-band data that we are working with is [Beeldmateriaal By using the `rast()` function along with the `lyrs` parameter, we can read specific raster bands; omitting this parameter would read instead all bands. -```r +``` r RGB_band1_TUD <- rast("data/tudlib-rgb.tif", lyrs = 1) ``` We need to convert this data to a data frame in order to plot it with `ggplot`. -```r +``` r RGB_band1_TUD_df <- as.data.frame(RGB_band1_TUD, xy = TRUE) ggplot() + @@ -68,13 +68,13 @@ This raster contains values between 0 and 255. These values represent degrees of We can use the `rast()` function to import specific bands in our raster object by specifying which band we want with `lyrs = N` (N represents the band number we want to work with). To import the green band, we would use `lyrs = 2`. -```r +``` r RGB_band2_TUD <- rast("data/tudlib-rgb.tif", lyrs = 2) ``` We can convert this data to a data frame and plot the same way we plotted the red band: -```r +``` r RGB_band2_TUD_df <- as.data.frame(RGB_band2_TUD, xy = TRUE) ggplot() + @@ -91,17 +91,17 @@ Next, we will work with all three image bands (red, green and blue) as an R rast To bring in all bands of a multi-band raster, we use the `rast()` function. -```r +``` r RGB_stack_TUD <- rast("data/tudlib-rgb.tif") ``` Let’s preview the attributes of our stack object: -```r +``` r RGB_stack_TUD ``` -```output +``` output class : SpatRaster dimensions : 4988, 4866, 3 (nrow, ncol, nlyr) resolution : 0.08, 0.08 (x, y) @@ -113,11 +113,11 @@ names : tudlib-rgb_1, tudlib-rgb_2, tudlib-rgb_3 ``` We can view the attributes of each band in the stack in a single output. For example, if we had hundreds of bands, we could specify which band we’d like to view attributes for using an index value: -```r +``` r RGB_stack_TUD[[2]] ``` -```output +``` output class : SpatRaster dimensions : 4988, 4866, 1 (nrow, ncol, nlyr) resolution : 0.08, 0.08 (x, y) @@ -129,17 +129,17 @@ name : tudlib-rgb_2 We can also use `ggplot2` to plot the data in any layer of our raster object. Remember, we need to convert to a data frame first. -```r +``` r RGB_stack_TUD_df <- as.data.frame(RGB_stack_TUD, xy = TRUE) ``` Each band in our RasterStack gets its own column in the data frame. Thus we have: -```r +``` r str(RGB_stack_TUD_df) ``` -```output +``` output 'data.frame': 24271608 obs. of 5 variables: $ x : num 85272 85272 85272 85272 85272 ... $ y : num 446694 446694 446694 446694 446694 ... @@ -149,7 +149,7 @@ str(RGB_stack_TUD_df) ``` Let’s create a histogram of the first band: -```r +``` r ggplot() + geom_histogram(data = RGB_stack_TUD_df, aes(`tudlib-rgb_1`)) ``` @@ -157,7 +157,7 @@ ggplot() + And a raster plot of the second band: -```r +``` r ggplot() + geom_raster(data = RGB_stack_TUD_df, aes(x = x, y = y, alpha = `tudlib-rgb_2`)) + @@ -180,7 +180,7 @@ This function allows us to: Let’s plot our 3-band image. Note that we can use the `plotRGB()` function directly with our RasterStack object (we don’t need a data frame as this function isn’t part of the `ggplot2` package). -```r +``` r plotRGB(RGB_stack_TUD, r = 1, g = 2, b = 3) ``` diff --git a/18-import-and-visualise-osm-data.md b/18-import-and-visualise-osm-data.md index 9480517b..bcc878d2 100644 --- a/18-import-and-visualise-osm-data.md +++ b/18-import-and-visualise-osm-data.md @@ -37,7 +37,7 @@ The geospatial data underlying this interface is made of geometrical objects (i. ## How to extract geospatial data from OpenStreetMap? -```r +``` r library(tidyverse) library(sf) assign("has_internet_via_proxy", TRUE, environment(curl::has_internet)) @@ -59,14 +59,14 @@ Beware that downloading and analysing the data for larger cities might be long, We first geocode our spatial text search and extract the corresponding bounding box (`getbb`). -```r +``` r library(osmdata) bb <- osmdata::getbb("Brielle") bb ``` -```output +``` output min max x 4.135294 4.229222 y 51.883500 51.931410 @@ -79,7 +79,7 @@ y 51.883500 51.931410 If you encounter an error linked to your internet proxy ("Error: Overpass query unavailable without internet R"), run this line of code. It might not be needed, but ensures that your machine knows it has internet. -```r +``` r assign("has_internet_via_proxy", TRUE, environment(curl::has_internet)) ``` @@ -98,12 +98,12 @@ For example: Brielle (Netherlands) and Brielle (New Jersey) By default, `getbb()` from the `osmdata` package returns the first item. This means that regardless of the number of returned locations with the given name, the function will return a bounding box and it might be that we are not looking for the first item. We should therefore try to be as unambiguous as possible by adding a country code or district name. -```r +``` r bb <- getbb("Brielle, NL") bb ``` -```output +``` output min max x 4.135294 4.229222 y 51.883500 51.931410 @@ -142,7 +142,7 @@ On this page we can read about the arguments needed for each function: a boundin -```r +``` r x <- opq(bbox = bb) %>% add_osm_feature(key = 'building') %>% osmdata_sf() @@ -157,11 +157,11 @@ What is this `x` object made of? It is a data frame of all the buildings contain -```r +``` r str(x$osm_polygons) ``` -```output +``` output Classes 'sf' and 'data.frame': 10632 obs. of 62 variables: $ osm_id : chr "40267783" "40267787" "40267791" "40267800" ... $ name : chr NA "F-51" "F-61" "F-7" ... @@ -248,7 +248,7 @@ Let's map the building age of post-1900 Brielle buildings. First, we are going to select the polygons and reproject them with the Amersfoort/RD New projection, suited for maps centred on the Netherlands. This code for this projection is: 28992. -```r +``` r buildings <- x$osm_polygons %>% st_transform(.,crs=28992) ``` @@ -268,7 +268,7 @@ Then we create a new variable using the threshold at 1900. Every date before 190 Then we use the `ggplot()` function to visualise the buildings by age. The specific function to represent information as a map is `geom_sf()`. The rest works like other graphs and visualisation, with `aes()` for the aesthetics. -```r +``` r start_date <- as.numeric(buildings$start_date) buildings$build_date <- if_else(start_date < 1900, 1900, start_date) @@ -291,7 +291,7 @@ We have produced a proof a concept on Brielle, but can we factorise our work to We might replace the name in the first line and run everything again. Or we can create a function. -```r +``` r extract_buildings <- function(cityname, year=1900){ nominatim_polygon <- geo_lite_sf(address = cityname, points_only = FALSE) bb <- st_bbox(nominatim_polygon) @@ -318,16 +318,16 @@ extract_buildings <- function(cityname, year=1900){ extract_buildings("Brielle, NL") ``` -```error +``` error Error in geo_lite_sf(address = cityname, points_only = FALSE): could not find function "geo_lite_sf" ``` -```r +``` r #test on Naarden extract_buildings("Naarden, NL") ``` -```error +``` error Error in geo_lite_sf(address = cityname, points_only = FALSE): could not find function "geo_lite_sf" ``` @@ -347,7 +347,7 @@ Leaflet is a ["open-source JavaScript library for mobile-friendly interactive ma ## One solution -```r +``` r #install.packages("leaflet") library(leaflet) @@ -370,8 +370,8 @@ leaflet(buildings2) %>% bringToFront = TRUE)) ``` -
- +
+ ::::::::::::::::::::::::::::::::: diff --git a/19-basic-gis-with-r-sf.md b/19-basic-gis-with-r-sf.md index bdbc3590..8fbf4be5 100644 --- a/19-basic-gis-with-r-sf.md +++ b/19-basic-gis-with-r-sf.md @@ -24,7 +24,7 @@ After completing this episode, participants should be able to… -```r +``` r library(tidyverse) library(sf) library(osmdata) @@ -48,7 +48,7 @@ Let's focus on old buildings and imagine we're in charge of their conservation. Let's select them and see where they are. -```r +``` r bb <- osmdata::getbb("Brielle, NL") x <- opq(bbox = bb) %>% add_osm_feature(key = 'building') %>% @@ -60,12 +60,12 @@ buildings <- x$osm_polygons %>% summary(buildings$start_date) ``` -```output +``` output Length Class Mode 10632 character character ``` -```r +``` r old <- 1800 # year prior to which you consider a building old buildings$start_date <- as.numeric(buildings$start_date) @@ -87,7 +87,7 @@ old_buildings <- buildings %>% If you encounter an error linked to your internet proxy ("Error: Overpass query unavailable without internet R"), run this line of code. It might not be needed, but ensures that your machine knows it has internet. -```r +``` r assign("has_internet_via_proxy", TRUE, environment(curl::has_internet)) ``` @@ -100,14 +100,14 @@ As conservationists, we want to create a zone around historical buildings where Let's say the conservation zone should be 100 meters. In GIS terms, we want to create a _buffer_ around polygons. The corresponding `sf` function is `st_buffer()`, with 2 arguments: the polygons around which to create buffers, and the radius of the buffer. -```r +``` r distance <- 100 # in meters #First, we check that the "old_buildings" layer projection is measured in meters: st_crs(old_buildings) ``` -```output +``` output Coordinate Reference System: User input: EPSG:28992 wkt: @@ -151,7 +151,7 @@ PROJCRS["Amersfoort / RD New", ID["EPSG",28992]] ``` -```r +``` r #then we use `st_buffer()` buffer_old_buildings <- st_buffer(x = old_buildings, dist = distance) @@ -168,7 +168,7 @@ ggplot(data = buffer_old_buildings) + Now, we have a lot of overlapping buffers. We would rather create a unique conservation zone rather than overlapping ones in that case. So we have to fuse the overlapping buffers into one polygon. This operation is called _union_ and the corresponding function is `st_union()`. -```r +``` r single_old_buffer <- st_union(buffer_old_buildings) %>% st_cast(to = "POLYGON") %>% st_as_sf() @@ -186,7 +186,7 @@ We create unique IDs to identify the new polygons. For the sake of visualisation speed, we would like to represent buildings by a single point (for instance: their geometric centre) rather than their actual footprint. This operation means defining their _centroid_ and the corresponding function is `st_centroid()`. -```r +``` r sf::sf_use_s2(FALSE) # s2 works with geographic projections, so to calculate centroids in projected CRS units (meters), we need to disable it. centroids_old <- st_centroid(old_buildings) %>% @@ -205,7 +205,7 @@ Now, we would like to distinguish conservation areas based on the number of hist -```r +``` r centroids_buffers <- st_intersection(centroids_old,single_old_buffer) %>% mutate(n = 1) @@ -239,7 +239,7 @@ We aggregate them by ID number (`group_by(ID)`) and sum the variable `n` to know Let's map this layer over the initial map of individual buildings, and save the result. -```r +``` r p <- ggplot() + geom_sf(data = buildings) + geom_sf(data = single_buffer, aes(fill=n_buildings), colour = NA) + @@ -255,7 +255,7 @@ p <- ggplot() + -```r +``` r ggsave(filename = "fig/ConservationBrielle.png", plot = p) ``` @@ -270,7 +270,7 @@ The historical threshold now applies to all pre-war buildings, but the distance :::::::::::::::::::::::: solution -```r +``` r old <- 1939 distance <- 10 @@ -320,7 +320,7 @@ pnew <- ggplot() + -```r +``` r ggsave(filename = "fig/ConservationBrielle_newrules.png", plot = pnew) ``` @@ -335,7 +335,7 @@ ggsave(filename = "fig/ConservationBrielle_newrules.png", ## Area -```r +``` r single_buffer$area <- sf::st_area(single_buffer) %>% units::set_units(., km^2) diff --git a/md5sum.txt b/md5sum.txt index e742e9d4..e3873793 100644 --- a/md5sum.txt +++ b/md5sum.txt @@ -4,7 +4,7 @@ "config.yaml" "eba2e90f9f89b711e953fb45172634ef" "site/built/config.yaml" "2024-07-02" "index.md" "db222f765de8e62e7990ec8ba4328199" "site/built/index.md" "2024-07-02" "links.md" "8184cf4149eafbf03ce8da8ff0778c14" "site/built/links.md" "2024-07-02" -"renv.lock" "e9277c61663340dc9c31ae9fe9347500" "site/built/renv.lock" "2024-07-02" +"renv.lock" "c307f5f1316bb5a3ac59513e62a124d7" "site/built/renv.lock" "2024-07-02" "episodes/01-intro-to-r.Rmd" "b13d8b2857351dc5b861e6f3f86f8add" "site/built/01-intro-to-r.md" "2024-07-02" "episodes/02-data-structures.Rmd" "2e3bbbc3ff746151b6bfeffdce1d61c9" "site/built/02-data-structures.md" "2024-07-02" "episodes/03-explore-data.Rmd" "9e557be2898dfa30ff101cb2717b7deb" "site/built/03-explore-data.md" "2024-07-02" @@ -25,4 +25,4 @@ "learners/reference.md" "1c7cc4e229304d9806a13f69ca1b8ba4" "site/built/reference.md" "2024-07-02" "learners/setup.md" "a3f579bba088d6c799abcf53dffb8aff" "site/built/setup.md" "2024-07-02" "profiles/learner-profiles.md" "60b93493cf1da06dfd63255d73854461" "site/built/learner-profiles.md" "2024-07-02" -"renv/profiles/lesson-requirements/renv.lock" "e9277c61663340dc9c31ae9fe9347500" "site/built/renv.lock" "2024-07-02" +"renv/profiles/lesson-requirements/renv.lock" "c307f5f1316bb5a3ac59513e62a124d7" "site/built/renv.lock" "2024-07-02"