mapnight1.qmd

# Lab 1 {#sec-mapnight1 .unnumbered}

We will start by exploring raster data processing and choropleths we can make from it. Raster data, consisting of gridded cells, allows us to represent continuous geographic phenomena such as temperature, elevation, or satellite imagery. Choropleths, on the other hand, are an effective way to visualize spatial patterns through the use of color-coded regions, making them invaluable for displaying discrete data like population density or election results. By combining these techniques, we will gain a comprehensive toolkit for conveying complex geographical information in a visually compelling manner.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# Provides various utility functions for R programming.
library(R.utils)
# For data manipulation and transformation.
library(dplyr)
# Spatial data
library(sf)
# Popular data visualization package in R.  
library(ggplot2)
# For creating thematic maps 
library(tmap)
# Color palettes suitable for data visualization, especially for those with color vision deficiencies.
library(viridis)
# A collection of color palettes for data visualization.
library(RColorBrewer)
# For working with raster data, such as gridded spatial data like satellite imagery or elevation data.
library(raster)
# An alternative to the 'raster' package and is used for working with large raster datasets efficiently.
library(terra)
# Tools for extracting data from raster layers at exact locations, often used in spatial analysis.
library(exactextractr)
# Common methods of the tidyverse packages for objects created with the {terra} package: SpatRaster and SpatVector
library(tidyterra)
# Querying Open Street Map data
library(osmdata)
```

## NOAA Night Lights data

### Import raster data

```{r import_a_raster}
MENA_lights <- rast("data/MENA_noaa_projected.tif")
```

Plot it.

```{r plot}
plot(MENA_lights) 
```

Have a look at the CRS.

```{r checkcrs}
crs(MENA_lights)
```

### Import the MENA shapefile

Import the Middle East and North Africa (MENA) shapefile, plot it, and verify its Coordinate Reference System (CRS). Is it the same as the raster's CRS?

```{r import_shp}
MENA_adm1 <- read_sf("data/MENA_projected.shp")
plot(MENA_adm1$geometry)
crs(MENA_adm1)
```

### Reproject the Raster

As we are using both the `raster` and `terra` packages to handle the raster data it is useful to write `terra::` or `raster::` in front of the function we are using.

We use the terra `project()` function, we need to define two things:

1.  The object we want to reproject and
2.  The CRS that we want to reproject it to.

```{r reproject}
MENA_lights <- terra::project(MENA_lights, crs(MENA_adm1)) # reporjectig the raster data to the crs of the MENA shapefile
crs(MENA_lights)
```

### Cropping and Masking

Cropping and masking are both spatial operations used in raster data analysis.

**Cropping**:

-   Purpose: Cropping a raster involves changing the extent of the raster dataset by specifying a new bounding box or geographic area of interest. The result is a new raster that covers only the specified region.

-   Typical Use: Cropping is commonly used when you want to reduce the size of a raster dataset to focus on a smaller geographic area of interest while retaining all the original data values within that area.

**Masking**:

-   Purpose: Applying a binary mask to the dataset. The mask is typically a separate raster or polygon layer where certain areas are designated as "masked" (1) or "unmasked" (0).

-   Typical Use: Masking is used when you want to extract or isolate specific areas or features within a raster dataset. For example, you might use a mask to extract land cover information within the boundaries of a protected national park.

In many cases, these cropping and masking are executed one after the other because it is computationally easier to crop when dealing with large datasets, and then masking.

```{r crop}
MENA_lights_crop <- crop(MENA_lights, extent(MENA_adm1))
MENA_lights_mask <- mask(MENA_lights_crop, MENA_adm1)
```

### Simple plot 

```{r plottogether}
plot(MENA_lights_mask)
plot(MENA_adm1$geometry, col= NA, add=T)
```

Let's improve this a bit. Remember that there is a lot we can do with [ColorBrewer](https://colorbrewer2.org/#type=sequential&scheme=BuGn&n=3).

::: {.panel-tabset group="data"}

# `plot`

```{r changecolours}
#| warning = FALSE

pal = rev(brewer.pal(6,"YlGnBu"))
plot(MENA_lights_mask, breaks=c(0,10,20,30,40,Inf), col=pal)
plot(MENA_adm1$geometry, col= NA, add=T)
```

# `tmap`

```{r}
#| warning = FALSE

# Define the palette
pal = rev(brewer.pal(6,"YlGnBu"))

# Create the base map
tm_shape(MENA_lights_mask) +
  tm_raster(breaks = c(0,10,20,30,40,Inf),
            palette = pal) +  # Plot the raster with breaks and palette
  tm_shape(MENA_adm1) +
  tm_borders(lwd = 1, col = "darkgrey") +  # Add borders to the administrative boundaries
tm_layout(frame = FALSE, legend.outside = TRUE, legend.outside.position = "right", 
          title = "Night light intensity")  # Remove frame
```
:::

::: {.callout-note title="Questions"}

Questions to ask yourself about how you can improve these maps, going back to [geo-visualisation and choropleths](https://pietrostefani.github.io/gds/mapvector.html).

-   What are the logical breaks for night lights data?

-   What should the colours be?

-   Have a look at some of the [`tmap` documentation](https://r-tmap.github.io/tmap/reference/tm_layout.html) to improve your map further .

:::

### Spatial join with vector data

Keep Egypt from all Middle East and North Africa Countries.

```{r egypt}
# Subset MENA_adm1 to select only Egypt
Egypt <- MENA_adm1 %>%
  filter(name == "Egypt")

# Crop MENA_adm1 raster to Egypt
Egypt_lights_mask <- crop(MENA_lights_mask, extent(Egypt))
```

Using `osmdata` package to query amenity data. We only query hospitals to keep it brief. 

```{r osm}
# creating bounding box for Cairo
egypt_bb <- getbb("Cairo Egypt")
cairo <- opq(bbox = egypt_bb)

# Getting hospitals in Cairo
amenities <- cairo %>%
   add_osm_feature(key = "amenity", value = c("hospital")) %>%
   osmdata_sf ()

# checking object and types of data
amenities
```

Keeping just the OSM points from the OSM data.

```{r just osm points}
amenities_points <- amenities$osm_points
```

## Plotting

An initial plot.
```{r plotting}
# Plotting
ggplot() + 
  geom_sf(data = amenities$osm_points, fill = 'orange') + 
  geom_spatraster(data=Egypt_lights_mask, alpha = 0.5) +
  coord_sf(xlim = c(30.852356, 31.717529), ylim = c(29.699982, 30.473532)) + theme_minimal()
```

Starting to improve the plot.
``` {r plottingmore}
#| warning = FALSE
# Plotting
ggplot() + 
  geom_sf(data = amenities$osm_points, aes(color = "Hospitals"), fill = 'red', size = 1, stroke = 1, color = "darkred") + 
  geom_spatraster(data = Egypt_lights_mask, aes(fill = ..value..), alpha = 0.5) +
  scale_fill_gradient(low = "black", high = "yellow", name = "Night Lights") +
  coord_sf(xlim = c(30.852356, 31.717529), ylim = c(29.699982, 30.473532)) + 
  labs(title = "Distribution of Hospitals and Night Lights in Egypt",
       x = "Longitude", y = "Latitude") +
  theme_minimal() +
  theme(legend.position = "bottom")

```

More on `geom_spatraster` (here)[https://dieghernan.github.io/tidyterra/reference/geom_spatraster.html].

## Extracting

You might want to extract values from a raster data set, and then map them in a vector `sf` framework or extract them to analyse them statistically. If it therefore very useful to know how to extract:

```{r extract}
#| warning = FALSE
# Using the 'raster::extract' function, it calculates the illumination values at the coordinates of the points. 
lights_hospitals <- raster::extract(Egypt_lights_mask,
                                amenities_points)

# Attach nightlight data at each point to the amenities dataframe
lights_hospitals <- cbind(amenities_points, lights_hospitals)

# Keep only specified columns
lights_hospitals <- lights_hospitals %>%
  select(MENA_noaa_projected, ID, geometry, osm_id)

# Check out the data
head(lights_hospitals)
```

::: callout-important
Make sure all your data is in the same CRS, otherwise the `raster::extract` function will not work properly.

You should not be seeing NAs, if you do you should use the terra `project()` function.

There is also the `exactextractor` package which we will use in the next section.
:::