forked from geocompx/geocompr
-
Notifications
You must be signed in to change notification settings - Fork 0
/
_02-ex.Rmd
115 lines (95 loc) · 4.64 KB
/
_02-ex.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
```{r 02-ex-e0, message=FALSE}
library(sf)
library(spData)
library(terra)
```
E1. Use `summary()` on the geometry column of the `world` data object that is included in the **spData** package. What does the output tell us about:
- Its geometry type?
- The number of countries?
- Its coordinate reference system (CRS)?
```{r 02-ex-e1}
summary(world)
# - Its geometry type?
# multipolygon
# - The number of countries?
# 177
# - Its coordinate reference system (CRS)?
# epsg:4326
```
E2. Run the code that 'generated' the map of the world in Section 2.2.3 (Basic map making).
Find two similarities and two differences between the image on your computer and that in the book.
- What does the `cex` argument do (see `?plot`)?
- Why was `cex` set to the `sqrt(world$pop) / 10000`?
- Bonus: experiment with different ways to visualize the global population.
```{r 02-ex-e2}
plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)
plot(st_geometry(world_cents), add = TRUE, cex = cex)
# - What does the `cex` argument do (see `?plot`)?
# It specifies the size of the circles
# - Why was `cex` set to the `sqrt(world$pop) / 10000`?
# So the circles would be visible for small countries but not too large for large countries, also because area increases as a linear function of the square route of the diameter defined by `cex`
# - Bonus: experiment with different ways to visualize the global population.
plot(st_geometry(world_cents), cex = world$pop / 1e9)
plot(st_geometry(world_cents), cex = world$pop / 1e8)
plot(world["pop"])
plot(world["pop"], logz = TRUE)
# Similarities: global extent, colorscheme, relative size of circles
#
# Differences: projection (Antarctica is much smaller for example), graticules, location of points in the countries.
#
# To understand these differences read-over, run, and experiment with different argument values in this script: https://github.com/Robinlovelace/geocompr/raw/main/code/02-contpop.R
#
# `cex` refers to the diameter of symbols plotted, as explained by the help page `?graphics::points`. It is an acronym for 'Chacter symbol EXpansion'.
# It was set to the square route of the population divided by 10,000 because a) otherwise the symbols would not fit on the map and b) to make circle area proportional to population.
```
E3. Use `plot()` to create maps of Nigeria in context (see Section 2.2.3).
- Adjust the `lwd`, `col` and `expandBB` arguments of `plot()`.
- Challenge: read the documentation of `text()` and annotate the map.
```{r 02-ex-e3}
nigeria = world[world$name_long == "Nigeria", ]
plot(st_geometry(nigeria), expandBB = c(0, 0.2, 0.1, 1), col = "gray", lwd = 3)
plot(world[0], add = TRUE)
world_coords = st_coordinates(world_cents)
text(world_coords, world$iso_a2)
# Alternative answer:
nigeria = world[world$name_long == "Nigeria", ]
africa = world[world$continent == "Africa", ]
plot(st_geometry(nigeria), col = "white", lwd = 3, main = "Nigeria in context", border = "lightgrey", expandBB = c(0.5, 0.2, 0.5, 0.2))
plot(st_geometry(world), lty = 3, add = TRUE, border = "grey")
plot(st_geometry(nigeria), col = "yellow", add = TRUE, border = "darkgrey")
a = africa[grepl("Niger", africa$name_long), ]
ncentre = st_centroid(a)
ncentre_num = st_coordinates(ncentre)
text(x = ncentre_num[, 1], y = ncentre_num[, 2], labels = a$name_long)
```
E4. Create an empty `SpatRaster` object called `my_raster` with 10 columns and 10 rows.
Assign random values between 0 and 10 to the new raster and plot it.
```{r 02-ex-e4, message = FALSE}
my_raster = rast(ncol = 10, nrow = 10,
vals = sample(0:10, size = 10 * 10, replace = TRUE))
plot(my_raster)
```
E5. Read-in the `raster/nlcd.tif` file from the **spDataLarge** package.
What kind of information can you get about the properties of this file?
```{r 02-ex-e5, message = FALSE}
nlcd = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
dim(nlcd) # dimensions
res(nlcd) # resolution
ext(nlcd) # extent
nlyr(nlcd) # number of layers
cat(crs(nlcd)) # CRS
```
E6. Check the CRS of the `raster/nlcd.tif` file from the **spDataLarge** package.
What kind of information you can learn from it?
```{r 02-ex-e6, message = FALSE}
cat(crs(nlcd))
```
```{asis, message = FALSE}
The WKT above describes a two-dimensional projected coordinate reference system.
It is based on the GRS 1980 ellipsoid with North American Datum 1983 and the Greenwich prime meridian.
It used the Transverse Mercator projection to transform from geographic to projected CRS (UTM zone 12N).
Its first axis is related to eastness, while the second one is related to northness, and both axes have units in meters.
The SRID of the above CRS is "EPSG:26912".
```