##  {background-image="cover.jpeg"}

### Outlines

-   What is spatial data and why should we use it?

-   Types of spatial data

-   Spatial data in R: how to store, load, and tidy spatial data

-   Visualizing spatial data

-   Exploring spatial data: Spatial authorization and spatial statistics to quantify relationships

------------------------------------------------------------------------

## **What is Spatial data and why should we use it?**

-   [**Spatial data, also known as geospatial data**]{style="color:red;"}, refers to information that identifies the geographic location and characteristics of natural or constructed features and boundaries on the Earth.

-   This data is often represented in terms of [**Cartesian coordinates (x,y)**]{style="color:red ;"} for two-dimensional maps, but may also include [**altitude (z)**]{style="color:red ;"} for a three-dimensional representation.

-   Spatial data can come in various forms including [**points, lines, and polygons**]{style="color:red ;"}.

-   The concept of spatial data is integral to a variety of applications that require an understanding of how different elements relate to each other within a geographical space.

-   This data can be collected through various means, including but not limited to, **satellite imagery**, **aerial photography**, and **ground-based surveys**.

## Spatial Data Formats

- Let us see a basic way to represent the spatial data. 

- But there is a variety of data formats to represent the data to suit different applications. 

- In most cases, spatial data formats are an extension of existing data formats.

    | **Type**          | **Non-Spatial Data**       | **Spatial Data**                    |
    |------------------|------------------------|------------------------------|
    | Text              | csv, json, xml             | csv, geojson, gml, kml              |
    | Binary/Compressed | pdf, xls, zip              | **[shapefile]{style="color:red;"}**, geopdf, geopackage       |
    | Images            | tiff, jpg, png             | geotiff, jpeg2000                   |
    | Databases         | SQLite, PostgreSQL, Oracle | Spatialite, PostGIS, Oracle Spatial |

## Spatial Data Types

Spatial Data can be broadly categorized into 2 types - Vector and Raster.

| **Type**                               | **Sub Types**                                                                                                            | **Examples**                                                                                                     |
|------------------|----------------------------|--------------------------|
| [**Vector Data**]{style="color:blue;"} | [**Points:**]{style="color:red;"} Represents features with a single coordinate pair [**(x, y)**]{style="color:purple;"}. | Locations of [**ATMs**]{style="color:green;"}, [**tree**]{style="color:green;"}                                  |
|                                        | [**Lines:**]{style="color:red;"} Represents linear features as ordered sequences of points.                              | [**Roads**]{style="color:green;"}, [**rivers**]{style="color:green;"}, [**utility lines**]{style="color:green;"} |
|                                        | [**Polygons:**]{style="color:red;"} Represents areas enclosed by closed loops of lines.                                  | [**Buildings**]{style="color:green;"}, [**lakes**]{style="color:green;"}, [**zones**]{style="color:green;"}      |
| [**Raster Data**]{style="color:blue;"} | **Grids:** Represents continuous data across a surface.                                                                  | [**Satellite images**]{style="color:green;"}, [**digital elevation models (DEMs)**]{style="color:green;"}        |
|                                        | **Pixels:** Smallest units in a raster dataset, each with a specific value.                                              | Values representing color, temperature,                                                                          |

------------------------------------------------------------------------

## Other data associated with Vector and Raster

-   [**Attribute data:**]{style="color:blue;"} Additional information describing the characteristics of spatial features.
-   For example in Vector data:
    -   **Schools:** have school name, number of students, education level.
    -   **Roads:**have road type, traffic volume, maintenance status.
    -   **Regions:** have population density, average income, and land use type.

**In Raster Data:** Each pixel can have multiple attributes, such as vegetation type or pollution levels.

## What are the components of spatial data

-   [**Geometry:**]{style="color:#ae01c7;"} Refers to the coordinates that define the shape of an object.
    -   [**Coordinates:**]{style="color:blue;"} These can be in different coordinate systems, such as geographic ([latitude]{style="color:purple;"}) and [longitude]{style="color:purple;"}) or projected systems.
    -   [**Types of Geometries:**]{style="color:blue;"} Common types include [points]{style="color:blue;"}, [lines]{style="color:blue;"}, and [polygons]{style="color:blue;"}.
    -   [**Dimension:**]{style="color:blue;"} Geometry can also have dimensional attributes like [2D]{style="color:purple;"} (x, y), [3D]{style="color:purple;"} (x, y, z), or even [4D]{style="color:purple;"} (x, y, z, time).
-   [**Topology:**]{style="color:#ae01c7;"} Refers to the spatial relationships between geometric entities.
    -   **Spatial Relationships**
        -   *Adjacency*: How entities are next to each other.
        -   *Connectivity*: How entities are connected.
        -   *Containment*: How one entity is contained within another.
    -   **Rules and Constraints:** Enforces rules such as:
        -   No overlapping polygons.
        -   Lines meeting at nodes.
        -   Ensuring polygon boundaries are closed loops.

## 

[**Attribute Data:**]{style="color:#ae01c7;"} provide detailed information about spatial features.

-   Example: A point representing a city may have attributes like population, elevation, and city name.

-   **Types of Attributes:** Attribute data is either,

    -   **Numeric** (population size),
    -   **Categorical** (land use type), or
    -   **Temporal** (date of a satellite image).

-   Attribute data is often stored in tables and linked to geometry through unique identifiers.

## How does the workflow of handling spatial data effectively

**Data Collection**: Data can be collected in either

-   **Satellite Imagery**: Captured by remote sensing satellites, these images can provide data on land cover, vegetation, weather patterns, and more.

-   **GPS Surveys**: Use Global Positioning System technology to collect precise location data for mapping and navigation purposes.

-   **Traditional Surveying Methods**: Involve measuring angles, distances, and elevations to map out areas accurately, often using tools like theodolites and total stations.

------------------------------------------------------------------------

**Data Storage**:

-   **`Shapefiles`**: A popular vector data format that stores geometric locations and associated attribute data.
    -   It consists of `at least three files` with extensions `.shp` (geometry), `.shx` (shape index), and `.dbf` (attribute data).
-   **`GeoJSON`**: A format for encoding a variety of geographic data structures using JavaScript Object Notation (JSON).
-   **`GeoTIFF`**: A raster format that includes spatial metadata embedded within the file, making it suitable for georeferenced raster data like satellite images.
-   **`Spatial Databases`**: PostgreSQL with PostGIS extensions enables the storage, querying, and manipulation of spatial data within a relational database system, supporting complex spatial operations and large datasets.

## Practice for Spatial Data in R

-   **Spatial Data**: Information about the location and shape of geographic features.
-   **Types**: Vector Data and Raster Data.

### Shapefile

-   **Shapefile**: A common format for representing vector data.
    -   **`.shp`**: Contains geometry data.
    -   **`.shx`**: Positional index of the geometry data.
    -   **`.dbf`**: Stores attributes for each shape.
-   **Additional Files**:
    -   **`.prj`**: Describes the projection (plain text).
    -   **`.sbn`** and **`.sbx`**: Spatial indices of the geometry data.
    -   **`.shp.xml`**: Contains spatial metadata in XML format.

------------------------------------------------------------------------

## The sf package for spatial vector data

-   **The sf package** stores geometric features (termed simple features, hence sf) in a data frame.

-   Geographic data is stored in the special geometry column.

-   The **sf** package is a modern alternative to the traditional **sp**, **rgeos**, and **rgdal** packages.

-   Supports a wide range of geometry types: points, lines, polygons, and their 'multi' versions.

-   Useful for handling geographic vector data.

### Understanding Shapefiles

-   **Shapefiles** are a simple non-topological format to store geometric location and attribute information for geographical features.

-   They are commonly used in spatial analyses as a type of map.

-   We can read a shapefile or a `sf` object with the [`st_read()`](https://r-spatial.github.io/sf/reference/st_read.html) function of **sf**.

-   For example, here we read the `Eth_Region_2013.shp` shapefile of **sf** which contains the regions of Ethiopia, with their name, and Region and country code.

## Reading a Shapefile with `sf`

```{r}
library(sf)
ethR_shape <- st_read("Ethiopia_All/Eth_Region_2013.shp", 
                      quiet = TRUE)
class(ethR_shape)
head(ethR_shape, 2)
```

------------------------------------------------------------------------

-   The `sf` object `ethR_shape` is a `data.frame` containing a collection with 
     - **[11 simple features (rows)] {style="color:#ae01c7;"}** and [6 attributes (columns)]{style="color:blue;"} [plus a list-column with the geometry of each feature] {style="color:red;"}.

-   A `sf` object contains the following objects of class `sf`, `sfc` and `sfg`:
    -   `sf` (simple feature): *each row* of the `data.frame` [is a single simple feature consisting of attributes and geometry]{style="color:#ae01c7;"}.
    -   `sfc` (simple feature geometry list-column): the `geometry` *column* of the `data.frame` [is a list-column of class `sfc` with the geometry of each simple feature]{style="color:blue;"}.
    -   `sfg` (simple feature geometry): *each of the rows* of the `sfc` list-column corresponds to the simple feature geometry (`sfg`) of a single simple feature.

-   The **sf** package stores geometric features in a data frame.

-   Geographic data is stored in the special geometry column.

------------------------------------------------------------------------

-   Using `dim()` and `names()` functions we can see how many rows and columns it contains, as well as the names of the columns.

```{r}
dim(ethR_shape)
names(ethR_shape)
```

-   Each row represents a geographic element, in this case, one of the 11 regions of Ethiopia.
-   Columns contain information about the region, including the name and administrative code.
-   The last column is the geometry column.

## Displaying the First Few Geometries

```{r}
head(ethR_shape$geometry)
```

Each region is captured as a multipolygon geometry.

## Plotting the Shapefile

```{r}
eth_reg <- ethR_shape %>%
  dplyr::select(region_name = REGIONNAME)
plot(eth_reg, main = "Region name")
```

-   The default behavior of plot with an sf object is to create a map for each column of data that describes the geography.

## Plotting with ggplot2

```{r}
library(ggplot2)
region_map <- ggplot(ethR_shape) +geom_sf(aes(fill=REGIONNAME))+
  scale_fill_discrete("Region") + coord_sf(datum = NA) +
  theme_bw() + labs(x = NULL, y = NULL) + 
  ggtitle("Map of regions in Ethiopia")
print(region_map)
```

-   ggplot2 is flexible once you learn the syntax.
-   The plot shows each region filled with different colors.

## Adding Labels to the Map

-   Adding labels to each region ensures they do not overlap.

```{r}
map_label <- ggplot() + geom_sf(data = ethR_shape) +
  ggrepel::geom_label_repel(data = ethR_shape,
    aes(label = REGIONNAME, geometry = geometry),
    stat = "sf_coordinates",min.segment.length = 0) + 
  labs(x = NULL, y = NULL)+ coord_sf(datum = NA) +  theme_bw()
print(map_label)
```

-   The [`st_geometry()`](https://r-spatial.github.io/sf/reference/st_geometry.html) function can be used to retrieve the simple feature geometry list-column (`sfc`).

------------------------------------------------------------------------

```{r}
# Geometries printed in abbreviated form 
st_geometry(ethR_shape) # View complete geometry 
st_geometry(ethR_shape)[[1]]
```

## **Creating a `sf` object**

-   We can use the [`st_sf()`](https://r-spatial.github.io/sf/reference/sf.html) function to create a `sf` object by providing two elements, namely, a `data.frame` with the attributes of each feature, and a simple feature geometry list-column `sfc` containing simple feature geometries `sfg`.

-   In more detail, we create simple feature geometries `sfg` and use the [`st_sfc()`](https://r-spatial.github.io/sf/reference/sfc.html) function to create a simple feature geometry list-column `sfc` with them.

-   Then, we use [`st_sf()`](https://r-spatial.github.io/sf/reference/sf.html) to put the `data.frame` with the attributes and the simple feature geometry list-column `sfc` together.

Simple feature geometries `sfg` objects can be, for example, of type `POINT` (single point), `MULTIPOINT` (set of points) or `POLYGON` (polygon), and can be created with [`st_point()`](https://r-spatial.github.io/sf/reference/st.html), [`st_multipoint()`](https://r-spatial.github.io/sf/reference/st.html) and [`st_polygon()`](https://r-spatial.github.io/sf/reference/st.html), respectively.

------------------------------------------------------------------------

-   Here, we create a `sf` object containing two single points, a set of points, and a polygon, with one attribute.

-   First, we create the simple feature geometry objects (`sfg`) of type `POINT`, `MULTIPOINT`, and `POLYGON`.

-   Then, we use [`st_sfc()`](https://r-spatial.github.io/sf/reference/sfc.html) to create a simple feature geometry list-column `sfc` with the `sfg` objects.

-   Finally, we use [`st_sf()`](https://r-spatial.github.io/sf/reference/sf.html) to put the `data.frame` with the attribute and the simple feature geometry list-column `sfc` together.

------------------------------------------------------------------------

```{r}
# Single point (point as a vector) 
p1_sfg <- st_point(c(2, 2)) 
p2_sfg <- st_point(c(2.5, 3))  # Set of points 
p <- rbind(c(6, 2), c(6.1, 2.6), c(6.8, 2.5), 
           c(6.2, 1.5), c(6.8, 1.8)) 
mp_sfg <- st_multipoint(p)  
```

------------------------------------------------------------------------

```{r}
#|ccolumn: true
p1 <- rbind(c(10, 0), c(11, 0), c(13, 2), c(12, 4), 
            c(11, 4), c(10, 0)) 
p2 <- rbind(c(11, 1), c(11, 2), c(12, 2), c(11, 1)) 
pol_sfg <- st_polygon(list(p1, p2))  # Create sf object 
p_sfc <- st_sfc(p1_sfg, p2_sfg, mp_sfg, pol_sfg) 
df <- data.frame(v1 = c("A", "B", "C", "D"))

# Plot single points, set of points and polygon
p_sf <- st_sf(df, geometry = p_sfc)  
library(ggplot2) 
ggplot(p_sf) + geom_sf(aes(col = v1), size = 3) + 
  theme_bw()
```

------------------------------------------------------------------------

## **`st_*()` functions**

Common functions to manipulate `sf` objects include the following:

-   [`st_read()`](https://r-spatial.github.io/sf/reference/st_read.html) reads a `sf` object,

-   [`st_write()`](https://r-spatial.github.io/sf/reference/st_write.html) writes a `sf` object,

-   [`st_crs()`](https://r-spatial.github.io/sf/reference/st_crs.html) gets or sets a new coordinate reference system (CRS),

-   [`st_transform()`](https://r-spatial.github.io/sf/reference/st_transform.html) transforms data to a new CRS,

-   [`st_intersection()`](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) intersects `sf` objects,

-   [`st_union()`](https://r-spatial.github.io/sf/reference/geos_combine.html) combines several `sf` objects into one,

-   [`st_simplify()`](https://r-spatial.github.io/sf/reference/geos_unary.html) simplifies a `sf` object,

-   [`st_coordinates()`](https://r-spatial.github.io/sf/reference/st_coordinates.html) retrieves coordinates of a `sf` object,

-   [`st_as_sf()`](https://r-spatial.github.io/sf/reference/st_as_sf.html) converts a foreign object to a `sf` object.

------------------------------------------------------------------------

-   We can delete some of the polygons by taking a subset of the rows of `map`.

```{r}
# Delete polygon 
map <- ethR_shape %>%
  st_filter(ethR_shape[ethR_shape$REG_P_CODE %in% c("3", "4"), ])

ggplot(map) + geom_sf(aes(fill = HRname))  
```

------------------------------------------------------------------------

::: {layout-ncol="2"}
::: left
```{r}
# Combine geometries 
ggplot(st_union(map, by_feature=FALSE) %>% st_sf())+geom_sf() 
```
:::

::: right
```{r}
# Simplify
ggplot(st_simplify(map, dTolerance =10000))+geom_sf()
```
:::
:::

`sf` object obtained by deleting some of its polygons (top), combining polygons (middle), and simplifying polygons (bottom).

## Transforming Point Data to an sf Object

-   The `st_as_sf()` function allows us to convert a foreign object to an `sf` object.

-   This can be specifying in the argument `coords` the name of the columns that contain the point coordinates.

-   Example: Converting a Data Frame to an sf Object

-   Here, we use `st_as_sf()` to turn a data frame containing coordinates `long` and `lat` and two variables `place` and `value` into an `sf` object.

-   Then, we use `st_crs()` to set the coordinate reference system given by the EPSG code 4326 to represent longitude and latitude coordinates.

------------------------------------------------------------------------

```{r}
library(sf)
library(mapview)
d <- data.frame(
  place = c("Ethiopia", "Kenya", "Somalia", "Sudan"),
  long = c(39.9559, 37.9062, 45.0792, 30.2176),
  lat = c(9.145, -1.286389, 2.046934, 19.6133),
  value = c(200, 150, 100, 300))
class(d)

dsf <- st_as_sf(d, coords = c("long", "lat"))
st_crs(dsf) <- 4326
class(dsf)
mapview(dsf)
```

------------------------------------------------------------------------

## Counting the Number of Points within Polygons

### Using `st_intersects()`

-   We can use the `st_intersects()` function of `sf` to count the number of points within the polygons of an `sf` object.
-   The returned object is a list with feature ids intersected in each of the polygons.
-   We can use the `lengths()` function to calculate the number of points inside each feature.

### Example: Counting Points within Polygons

In this example, we create a map with divisions (an `sf` object) and generate random points over the map.

------------------------------------------------------------------------

We then count the number of points within each polygon using `st_intersects()` and visualize the results with `ggplot2`.

```{r}
library(sf)
library(ggplot2)
# Points over map (simple feature geometry list-column sfc)
points <- st_sample(map, size = 100)
# Map of points within polygons
ggplot() + geom_sf(data = map) + geom_sf(data = points)
# Intersection (first argument map, then points)
inter <- st_intersects(map, points)
# Add point count to each polygon
map$count <- lengths(inter)
# Map of number of points within polygons
ggplot(map) + geom_sf(aes(fill = count))
```

------------------------------------------------------------------------

## Identifying Polygons Containing Points

### Using `st_intersects()`

-   Given an `sf` object with points and an `sf` object with polygons, we can use the `st_intersects()` function to obtain the polygon each of the points belongs to.

### Example: Identifying Polygons Containing Points

-   In this example, we create a map with divisions (an `sf` object) and generate three random points over the map.

-   We then identify which polygons contain these points and add the polygon names to the points data.

```{r}
library(sf)
library(ggplot2)

# Points over map (sf object)
points <- st_sample(map, size = 3) %>% st_as_sf()
inter <- st_intersects(points, map)# first points, then map
points$areaname <- map[unlist(inter), "HRname", drop = TRUE] # drop geometry
points
# Map
ggplot(map) + geom_sf() + geom_sf(data = points) + 
 geom_sf_label(data = map[unlist(inter), ], aes(label = HRname), nudge_y = 0.2)
```

## Joining Map and Data

### Using `left_join()`

-   Sometimes, a map and its corresponding data are available separately and we may wish to create an `sf` object representing the map with the added data that we can manipulate and plot.

-   We can create an `sf` map with the data attributes by joining the map and the data with the `left_join()` function of the `dplyr` package.

### Example: Adding Air Pollution Data to a World Map

First, we use the `ne_countries()` function of `rnaturalearth` to download the world map with the country polygons of class `sf`.

```{r}
library(rnaturalearth)
map <- ne_countries(returnclass = "sf")

```

------------------------------------------------------------------------

-   Then, we use the `wbstats` package to download a data frame of air pollution data from the World Bank.
    -   Specifically, we search the pollution indicators with `wb_search()`, and use `wb_data()` to download PM2.5 in the year 2016 by specifying the indicator corresponding to PM2.5, and the start and end dates.

```{r}
library(wbstats)
indicators <- wb_search(pattern = "pollution")
#d <- wb_data(indicator = "EN.ATM.PM25.MC.M3", start_date = 2016, end_date = 2016)

```

------------------------------------------------------------------------

-   Next, we use the `left_join()` function of `dplyr` to join the map and the data, specifying the argument `by` with the variables we wish to join by. Here, we use the ISO3 standard code of the countries rather than the country names, since names can be written differently in the map and the data frame.

```{r, eval=FALSE}
library(dplyr)
library(ggplot2)
library(viridis)

map1 <- left_join(map, d, by = c("iso_a3" = "iso3c"))
ggplot(map1) + geom_sf(aes(fill = EN.ATM.PM25.MC.M3)) +
  scale_fill_viridis() + labs(fill = "PM2.5") + theme_bw()

```

------------------------------------------------------------------------

### Note

When we use `left_join()`, the class of the resulting object is the same as the class of the first argument.

```{r, eval=FALSE}
map1 <- left_join(map, d, by = c("iso_a3" = "iso3c"))
class(map1)
# [1] "sf" "data.frame"

d1 <- left_join(d, map, by = c("iso3c" = "iso_a3"))
class(d1)
# [1] "data.frame"

```

## Raster Data

-   **Raster Data**: A spatial data structure that divides the study region into equal-sized rectangles called cells or pixels, storing one or more values for each cell.
-   **Uses**: Represent spatially continuous phenomena such as elevation, temperature, or air pollution values.
-   **Main Packages**:
    -   **`terra`**: Primary package for working with raster data. Also supports vector data.
    -   **`raster`**: Previously used for raster data; **`terra`** is faster and has more functionality.
    -   **`stars`**: Used for analyzing raster data and spatial data cubes (arrays with spatial dimensions).

------------------------------------------------------------------------

## GeoTIFF

-   **GeoTIFF**: A common format for raster data, with the extension **`.tif`**.
-   **Example**: Reading the **`elev.tif`** file from the **`terra`** package, representing elevation in Luxembourg.

```{r}
# Load the terra package 
library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- terra::rast(pathraster)
print(r)
```

------------------------------------------------------------------------

```{r}
# Plot the raster data
plot(r)
```

------------------------------------------------------------------------

### EPSG Codes and CRS Transformation

#### EPSG Codes

-   Most common CRSs can be specified by providing their EPSG (European Petroleum Survey Group) codes or their Proj4 strings.

Details of a given projection can be inspected using the `st_crs()` function of the `sf` package.

-   Example: EPSG Code 4326

-   EPSG code 4326 refers to the WGS84 longitude/latitude projection.

```{r}
st_crs("EPSG:4326")$Name
st_crs("EPSG:4326")$proj4string
st_crs("EPSG:4326")$epsg
```

------------------------------------------------------------------------

## Transforming CRS with `sf` and `terra`

-   Functions `sf::st_crs()` and `terra::crs()` allow us to get the CRS of spatial data.

-   These functions also allow us to set a CRS to spatial data by using `st_crs(x) <- value` if `x` is a `sf` object, and `crs(r) <- value` if `r` is a raster.

### Important Note

-   Setting a CRS does not transform the data; it just changes the CRS label.

-   We may want to set a CRS to data that does not come with CRS, and the CRS should be what it is, not what we would like it to be.

-   We use `sf::st_transform()` and `terra::project()` to transform the `sf` or raster data, respectively, to a new CRS.

------------------------------------------------------------------------

## Working with `sf` Package

### Reading and Getting the CRS

```{r}
library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
map <- st_read(pathshp, quiet = TRUE)
# Get CRS
st_crs(map)
```

------------------------------------------------------------------------

## Transforming the CRS

```{r}
# Transform CRS
map2 <- st_transform(map, crs = "EPSG:4326")

# Get CRS of transformed data
st_crs(map2)

```

------------------------------------------------------------------------

## Working with `terra` Package

### Reading and Getting the CRS of a Raster

```{r}
library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
# Get CRS
crs(r)
```

------------------------------------------------------------------------

## Transforming the CRS of a Raster

```{r}
# Transform CRS
r2 <- terra::project(r, "EPSG:2169")
# Get CRS of transformed data
crs(r2)
```

------------------------------------------------------------------------

## The `terra` Package for Raster and Vector Data

The `terra` package (Hijmans 2022) provides functions to create, read, manipulate, and write raster and vector data. Raster data represents spatially continuous phenomena through a grid of equally sized cells, while vector data includes points, lines, and polygons with associated attributes. This chapter demonstrates how to handle raster and vector data using `terra`.

------------------------------------------------------------------------

## Raster Data

### Creating and Reading Raster Data

The `rast()` function can create and read raster data. The `writeRaster()` function allows writing raster data. Here, we read elevation data for Luxembourg from a file provided by `terra`.

```{r}
library(terra) 
pathraster <- system.file("ex/elev.tif", package = "terra") 
r <- rast(pathraster)
plot(r)
```

------------------------------------------------------------------------

We can also create a `SpatRaster` object by specifying dimensions and extents.

### Raster Operations

Several functions provide information about raster size and dimensions.

```{r}
nrow(r) # number of rows
ncol(r) # number of columns
dim(r) # dimensions
ncell(r) # number of cells
```

The `values()` function sets and accesses raster values.

```{r}
values(r) <- 1:ncell(r)

```

Creating multilayer rasters and subsetting layers is straightforward.

```{r}
r2 <- r * r
s <- c(r, r2)
plot(s[[2]]) # layer 2

```

------------------------------------------------------------------------

Generic operations on rasters include:

```{r}
plot(min(s))
plot(r + r + 10)
plot(round(r))
plot(r == 1)

```

------------------------------------------------------------------------

## Vector Data

The `SpatVector` class handles vector data with attributes. The `vect()` function reads shapefiles, and `writeVector()` writes `SpatVector` objects. Here, we obtain a map of Luxembourg's divisions.

```{r}
pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)

```

------------------------------------------------------------------------

Creating a `SpatVector` with point locations and attributes:

```{r}
long <- c(-0.118092, 2.349014, -3.703339, 12.496366)
lat <- c(51.509865, 48.864716, 40.416729, 41.902782)
longlat <- cbind(long, lat)

crspoints <- "+proj=longlat +datum=WGS84"
d <- data.frame(place = c("London", "Paris", "Madrid", "Rome"), value = c(200, 300, 400, 500))
pts <- vect(longlat, atts = d, crs = crspoints)
plot(pts)

```

------------------------------------------------------------------------

## Cropping, Masking, and Aggregating Raster Data

The `terra` package provides functions to crop, mask, and aggregate raster data. Here, we demonstrate these operations with temperature data for Spain.

### Downloading Data

```{r, eval=FALSE}
library(terra)
r <- geodata::worldclim_country(country = "Spain", var = "tavg", res = 10, path = tempdir())
plot(r)
```

------------------------------------------------------------------------

## Averaging Temperature Data

```{r, eval=FALSE}
r <- mean(r)
plot(r)

```

------------------------------------------------------------------------

## Cropping and Masking Data

```{r}
library(ggplot2)
library(terra)
map <- rnaturalearth::ne_states("Spain", returnclass = "sf")
map <- map[-which(map$region == "Canary Is."), ] # delete region
ggplot(map) + geom_sf()

```

------------------------------------------------------------------------

## Aggregating Data

```{r}
r <- terra::aggregate(r, fact = 20, fun = "mean", na.rm = TRUE)
plot(r)

```

------------------------------------------------------------------------

## Extracting Raster Values at Points

The `extract()` function retrieves raster values at specified points.

### Example

```{r}
library(terra)
r <- rast(system.file("ex/elev.tif", package = "terra"))
v <- vect(system.file("ex/lux.shp", package = "terra"))

points <- crds(centroids(v))

plot(r)
plot(v, add = TRUE)
points(points)

```

------------------------------------------------------------------------

## Extracting Values

```{r}
points <- as.data.frame(points)
valuesatpoints <- extract(r, points)
cbind(points, valuesatpoints)

```

## Extracting and Averaging Raster Values

Extracting raster values within polygons and computing area-weighted averages.

### Extracting Values Within Polygons

```{r}
# Extracted raster cells within each polygon
head(extract(r, v, na.rm = TRUE))
# Extracted raster cells and percentage of area
# covered within each polygon
head(extract(r, v, na.rm = TRUE, weights = TRUE))

```

------------------------------------------------------------------------

## Averaging Values

```{r}
# Average raster values by polygon
v$avg <- extract(r, v, mean, na.rm = TRUE)$elevation

# Area-weighted average raster values by polygon (weights = TRUE)
v$weightedavg <- extract(r, v, mean, na.rm = TRUE, weights = TRUE)$elevation

library(ggplot2)
library(tidyterra)

# Plot average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = avg)) + scale_fill_terrain_c()

# Plot area-weighted average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = weightedavg)) + scale_fill_terrain_c()

```

------------------------------------------------------------------------

## Introduction to Mapping in R

### Importance of Maps

-   Maps are essential tools for visualizing spatial data.
-   They help to reveal patterns and trends that are not immediately apparent in tabular data.
-   Useful in various fields such as geography, urban planning, epidemiology, and environmental science.

### Packages Covered

-   `ggplot2`: A powerful package for creating static maps based on the grammar of graphics.
-   `leaflet`: A package for creating interactive maps using the Leaflet JavaScript library.
-   `mapview`: A package for quick and easy interactive map visualization.
-   `tmap`: A versatile package that can create both static and interactive maps.
-   `flowmapblue`: A package for visualizing mobility flows with interactive maps.

------------------------------------------------------------------------

## Preparing Areal Data

### Data Source

-   The data used is from a study on sudden infant deaths in North Carolina, USA, for the years 1974 and 1979.

### Loading Data

-   Use the `sf` package to load shapefile data.

```{r}
library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
d$vble <- d$SID74
d$vble2 <- d$SID79
```

------------------------------------------------------------------------

## Static Maps with ggplot2

### Grammar of Graphics

-   `ggplot2` is based on the "grammar of graphics", which provides a structured way to describe and build graphs.

-   You define aesthetics (like colors and shapes) and layers (like points, lines, and polygons).

### Creating a Map

-   First, load the required libraries.

-   Use `geom_sf()` to plot the spatial data.

-   Customize the map with color scales and themes.

```{r}
library(ggplot2)
library(viridis)
ggplot(d) + geom_sf(aes(fill = vble)) +
  scale_fill_viridis() + theme_bw()
```

------------------------------------------------------------------------

### Saving Plots

-   Save your plots to a file using `ggsave()`.

    ## Interactive Maps with leaflet

    ### Creating a Leaflet Map

    -   Leaflet maps are interactive and require transforming the coordinate reference system (CRS) to EPSG:4326.

```{r}
d <- st_transform(d, 4326)

```

### Define Color Palette and Create Map

-   Define a color palette using `colorNumeric()`.

-   Create the map and add polygons with `addPolygons()`.

-   Add a legend with `addLegend()`.

    ```{r}
    library(leaflet)
    pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)
    leaflet(d) %>% addTiles() %>%
      addPolygons(color = "white", fillColor = ~ pal(vble),
                  fillOpacity = 0.8) %>%
      addLegend(pal = pal, values = ~vble, opacity = 0.8)

    ```

------------------------------------------------------------------------

## Quick Maps with mapview

### Creating a Map

-   Use `mapview()` for quick and simple interactive maps.

```{r}
library(mapview)
mapview(d, zcol = "vble")
```

### Customizing Maps

-   Customize the map by changing the background, color palette, and legend title.

```{r, eval=FALSE}
library(RColorBrewer)
palette <- colorRampPalette(brewer.pal(9, "YlOrRd"))

mapview(d, zcol = "vble", map.types = "CartoDB.DarkMatter",
        col.regions = colorRampPalette(brewer.pal(9, "YlOrRd"))(), 
        layer.name = "SDI")
```

## Side-by-Side and Synchronized Maps

### Side-by-Side Maps

-   Display multiple maps side-by-side for comparison.

```{r, eval=FALSE}
library(leaflet.extras2)
#m1 | m2
```

------------------------------------------------------------------------

### Synchronized Maps

-   Sync multiple maps to pan and zoom together using `leafsync`.

```{r, eval=FALSE}
m <- leafsync::sync(m1, m2)
```

## Static and Interactive Maps with tmap

### Static Map

-   Create static maps with `tmap` using `tm_shape()` and `tm_polygons()`.

```{r}
library(tmap)
tmap_mode("plot")
tm_shape(d) + tm_polygons("vble")
```

------------------------------------------------------------------------

### Interactive Map

-   Switch to interactive mode with `tmap_mode("view")` and create maps.

```{r}
tmap_mode("view")
tm_shape(d) + tm_polygons("vble")
```

------------------------------------------------------------------------

## Mapping Point Data

### Creating Point Data Map with ggplot2

-   Plot point data with `ggplot2`, using `geom_sf()` to define aesthetics for color and size.

```{r, eval=FALSE}
ggplot(d) + geom_sf(aes(col = vble, size = size)) +
  scale_color_viridis()
```

------------------------------------------------------------------------

## Mapping Raster Data

### Creating Raster Data Map with ggplot2

-   Plot raster data by converting it to a data frame and using `geom_raster()`.

```{r}
ggplot(d) + geom_sf() +
  geom_raster(data = as.data.frame(r, xy = TRUE),
    aes(x = x, y = y, fill = elevation))
```

------------------------------------------------------------------------

## Mapping Mobility Flows with flowmapblue

### Creating Interactive Flow Map

-   Create interactive flow maps with `flowmapblue`.

-   Define locations and flows data frames.

-   Generate the map with clustering, dark mode, and animation options.

```{r, eval=FALSE}
library(flowmapblue)
locations <- data.frame(
  id = c(1, 2, 3),
  name = c("New York", "London", "Rio de Janeiro"),
  lat = c(40.713543, 51.507425, -22.906241),
  lon = c(-74.011219, -0.127738, -43.180244))
flows <- data.frame(
  origin = c(1, 2, 3, 2, 1, 3),
  dest = c(2, 1, 1, 3, 3 , 2),
  count = c(42, 51, 50, 40, 22, 42))
flowmapblue(locations, flows, mapboxAccessToken,
            clustering = TRUE, darkMode = TRUE, animation = FALSE)
```

------------------------------------------------------------------------

## R Packages to Download Open Spatial Data

Spatial data are used in a wide range of disciplines including environment, health, agriculture, economy, and society (Moraga and Baker 2022). Several R packages have been recently developed as clients for various databases that can be used for easy access to spatial data including administrative boundaries, climatic, and OpenStreetMap data. Here, we give short reproducible examples on how to download and visualize spatial data that can be useful in different settings. More extended examples and details about the capabilities of each of the packages can be seen at the packages’ websites and the rspatialdata website which provides a collection of tutorials on R packages to download and visualize spatial data using R.

## Administrative Boundaries of Countries

### Overview

-   Administrative boundaries data can be obtained from multiple packages:
    -   `rnaturalearth`: Global administrative boundaries.
    -   `tidycensus` and `tigris`: USA-specific data.
    -   `mapSpain`: Spain-specific data.
    -   `geobr`: Brazil-specific data.
    -   `giscoR`: Eurostat GISCO data.

------------------------------------------------------------------------
