Point-in-polygon with R: Animated fire map of Africa
June 05, 2022 | Milos Popovic
Every summer we read sinister news about the spread of wildfires in Australia, Americas and Europe. For example, I still remember those terrible photos of the 2007 fires in Greece, which burned nearly as twice as the area of the 2021 fires and killed 84 people.
While fires are prevalent in the summer months and affect other areas of the globe as well, we hardly encounter any reports on African fires. And, African fires consume an area as large as Europe each year, emitting greenhouse gases and aerosols 😵.
In this tutorial, we’ll use the Near Real-Time (NRT) active fire data within 3 hours of satellite observation from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the satellites distributed from NASA’s Fire Information for Resource Management System (FIRMS). Then I’ll show you how to segment the national map of Africa into hex bins and identify which fire points fall within each of them. Finally, we’ll wrap up this tutorial with a timelapse map of African fires for the last 10 days.
Import libraries
In this tutorial, we use a few important libraries such as tidyverse and sf for spatial analysis and data wrangling; giscoR for importing the shapefile of Africa; data.table for downloading the data from an URL; as well as classInt for creating bins. We will also use gganimate to make a timelapse map of African active fires just as we did it here.
# libraries we need
libs <- c(
"tidyverse", "sf", "giscoR",
"data.table", "classInt",
"gganimate", "gifski",
"transformr"
)
# install missing libraries
installed_libs <- libs %in% rownames(installed.packages())
if (any(installed_libs == F)) {
install.packages(libs[!installed_libs])
}
# load libraries
invisible(lapply(libs, library, character.only = T))Fetch the Africa shapefile
We first import the national map of Africa in sf format using the usual suspect giscoR. We also store the bounding box of our Africa sf object as a vector because we will use this information to query the NASA FIRMS API later on.
# 1. GET THE AFRICA SHP
#----------------------
get_africa_sf <- function(africa) {
africa <- giscoR::gisco_get_countries(
year = "2020", epsg = "4326",
resolution = "10", region = "Africa"
)
return(africa)
}
africa <- get_africa_sf()
africa_coords <- st_bbox(africa) %>%
as.matrix() %>%
as.vector()Create a hex grid map
Since African countries and regions differ in size and larger polygons usually steal our attention due to their size, we will use a hexbin map to equalize the importance of each part of the map. We choose 100,000 m² for the size of the grid cell to better capture the spatial variation.
Given that the default degree-based coordinate reference system (CRS) of our Africa sf object doesn’t allow us to create a meaningful grid we first transform it into a metric-based CRS, i.e. 3575. Then we devise the formula to calculate the area size of each hexagon cell.
We intersect the grid and our Africa sf object and assign unique identifeir to each cell in order to be able to identify the area of interest. In the final step, we transform the object into multipolygons and embrace the degree-based CRS (EPSG:4326).
# 2. MAKE A HEXBIN MAP OF AFRICA
#-------------------------------
crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs"
get_africa_hex <- function(africa_transfomed, africa_hex, africa_final) {
africa_transfomed <- africa %>%
st_transform(3575)
# Make grid of 100,000 m2
africa_hex <- st_make_grid(africa_transfomed,
cellsize = (3 * sqrt(3) * 196.18875^2) / 2,
what = "polygons",
square = F
) %>%
st_intersection(africa_transfomed) %>%
st_sf() %>%
mutate(id = row_number()) %>%
filter(
st_geometry_type(.) %in% c("POLYGON", "MULTIPOLYGON")
) %>%
st_cast("MULTIPOLYGON") # transform to multipolygons
africa_final <- st_transform(africa_hex, crs = crsLONGLAT) # transform back to longlat
return(africa_final)
}
africa_final <- get_africa_hex()Retrieve NASA fire data
In this tutorial, we will access the NRT active fire data within 3 hours of satellite observation from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Aqua and Terra satellites, generously provided by NASA via this web service.
Once you land on the page of the web service you’ll notice a submit form with blank fields for the area of interest, satellite source, nominal date, day range and map key. In order to write a query we first need to retreive the map key. Scroll down until you reach the Map Key window. Then click on GET MAP_KEY button, enter your email address in the pop-up window and press GET MAP_KEY button once more. This will automatically send a map key to your designated email address.
We’ll need map key and several other components to build an URL from which we will be able to retrieve the active fires data. In the chunk below, we define a function that will do this for us.
# 3. GET THE FIRE DATA
#---------------------
get_fire_data <- function(main_url, map_key, source,
area_coords, day_range, date) {
url <- paste(main_url, map_key, source, area_coords,
day_range, date,
sep = "/"
)
fire_data <- data.table::fread(url)
return(fire_data)
}Now let’s define the parameters for the URL.
- The main_url will get us the data in CSV format;
- In map_key field you should replace the asterisks with the map key that you received via your email;
- The source indicates one of the three satellite sources and we choose the default option, MODIS NRT;
- Then there is area_coords where we call the bounding box of Africa in vector format and insert comma after every longitude/latitude value;
- Next, date_range is the number of days in the past from the nominal date that we would like to retrieve. In our case, we take the last 10 days, which is the maximum option;
- Finally, date is the nominal date and we take the yesterday’s date as the starting point in order to retrive the complete data. We plug all the components and call the function.
Enjoy the magic 🤩!
main_url <- "https://firms.modaps.eosdis.nasa.gov/api/area/csv"
map_key <- "********************************"
source <- "MODIS_NRT"
area_coords <- paste(africa_coords, sep = "", collapse = ",")
day_range <- 10
date <- Sys.Date() - 1
fire_data <- get_fire_data(
main_url = main_url, map_key = map_key, source = source,
area_coords = area_coords, day_range = day_range, date = date
)Below we access the fire data that we had just retrieved. Notice that apart from latitude, longitude and date (acq_date) we also retrieved several other fields. One of them, confidence, indicates the quality of individual hotspot/fire pixels in the range from 0 (no confidence) to 100% (complete confidence). In the following section, we’ll use this information to weigh every fire event.
head(fire_data)
latitude longitude brightness scan track acq_date acq_time satellite
31.64720 6.00213 310.59 1.17 1.08 2022-05-26 122 Aqua
31.64870 5.98990 305.19 1.17 1.08 2022-05-26 122 Aqua
31.79013 6.06903 305.27 1.17 1.07 2022-05-26 122 Aqua
32.94583 3.23200 303.61 1.74 1.29 2022-05-26 122 Aqua
32.95462 3.23890 303.11 1.74 1.29 2022-05-26 122 Aqua
27.76350 9.19572 303.12 1.04 1.02 2022-05-26 124 Aqua
instrument confidence version bright_t31 frp daynight
1: MODIS 80 6.1NRT 294.94 11.40 N
2: MODIS 63 6.1NRT 294.55 7.06 N
3: MODIS 64 6.1NRT 294.85 7.30 N
4: MODIS 45 6.1NRT 269.50 20.18 N
5: MODIS 35 6.1NRT 270.01 19.40 N
6: MODIS 53 6.1NRT 290.58 6.15 NUse point-in-polygon to filter African fires
When we applied the bounding box of Africa to our query we also captured fire events in the Arabian peninsula and South Europe. In the code below, we will use point-in-polygon technique to remove those active fires.
We transform our data.table object into a point sf object. You probably noticed that we call column 2 and then 1; that’s because longitude comes second and latitude comes first in our fire_data table! We also set EPSG:4326 CRS to match the one from the Africa sf object.
Then we apply function st_join to keep only those fire points that fall within our Africa sf object.
# 4. FILTER AND AGGREGATE AFRICA'S FIRE EVENTS
#---------------------------------------------
get_africa_fires <- function(fire_coords, fire_africa) {
fire_coords <- st_as_sf(fire_data, coords = c(2, 1)) %>%
st_set_crs(4326) %>%
st_transform(crsLONGLAT)
fire_africa <- st_join(fire_coords, africa_final, join = st_within)
return(fire_africa)
}
fire_africa <- get_africa_fires()Now we just have to aggregate the fire events on the cell/date levels. In the chunk below, we first weigh every fire event using the confidence column.
Then we turn the fire data into a data.frame and aggregate the weighted values on the cell/date levels. Next, we join the aggregated values on our hex sf object. In the last remaining steps, we clean up our data:
- we get rid of the decimal places by rounding the fire events;
- we turn the remaining 0s into missing values so as to omit them from our animated map
- and, finally, we declare acq_date to be date column so that gganimate is clear about the animation field.
get_aggregated_africa_fires <- function() {
fire_africa$weighted_fire <- fire_africa$confidence / 100
fire_africa_sum <- st_drop_geometry(fire_africa) %>%
group_by(id, acq_date) %>%
summarise_at(
vars(weighted_fire),
list(sum_fire = sum)
)
fire_africa_sf <- left_join(africa_final, fire_africa_sum, by = "id")
fire_africa_sf$sum_fire <- round(fire_africa_sf$sum_fire, 0)
fire_africa_sf$sum_fire[fire_africa_sf$sum_fire == 0] <- NA
fire_africa_sf$acq_date <- as.Date(fire_africa_sf$acq_date)
return(fire_africa_sf)
}
fire_africa_sf <- get_aggregated_africa_fires()Mapping Africa’s fires
We’re almost ready to map Africa’s active fires! Before doing so, we create a discrete scale for our measure because it’s much easier to distinguish among the colors when you have clear-cut categories and labels.
We first compute the natural interval of our measure to get our cutpoints. I usually stick to 4-6 cutpoints to keep the colors discernible. Then we create labels based on the cutpoints and round them to get rid of the decimals. Finally, we carve out the categorical variable cat based on the cutpoints.
# 5. ANIMATE
#-----------
get_intervals <- function(ni, labels, df) {
df <- drop_na(fire_africa_sf)
ni <- classIntervals(df$sum_fire,
n = 5,
style = 'jenks')$brks
# create categories
labels <- c()
for(i in 1:length(ni)){
labels <- c(labels, paste0(round(ni[i], 0),
"–",
round(ni[i + 1], 0)))
}
labels <- labels[1:length(labels)-1]
df$cat <- cut(df$sum_fire,
breaks = ni,
labels = labels,
include.lowest = T)
return(df)
}
df <- get_intervals()Let’s make that map, shall we?
In the following chunk, we cleanse our final data from any remaining missing values and plug it into ggplot. We fill the values using the carved categories and interact them with the observed date to make sure that fire events do not “move around” as the dates shift.
We will also use an overlay map of African nations that we had created at the outset. The fire events will be filled using the manually generated color values that resemble a well-known Plasma palette.
get_africa_map <- function(africa_map) {
africa_map <- ggplot(na.omit(df)) +
geom_sf(
mapping = aes(
fill = cat,
group = interaction(cat, acq_date)
),
color = NA, size = 0
) +
geom_sf(
data = africa,
fill = NA,
color = "white",
size = 0.05
) +
scale_fill_manual(
name = "",
values = rev(c(
"#3d0965", "#8a335c", "#c26958",
"#e3aa61", "#f1ef75"
)),
drop = F
) +
coord_sf(crs = crsLONGLAT) +
guides(fill = guide_legend(
direction = "horizontal",
keyheight = unit(1.5, units = "mm"),
keywidth = unit(12.5, units = "mm"),
title.position = "top",
title.hjust = .5,
label.hjust = .5,
nrow = 1,
byrow = T,
reverse = F,
label.position = "bottom"
)) +
theme_minimal() +
theme(
axis.line = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.position = c(.5, .125),
legend.text = element_text(size = 20, color = "white"),
panel.grid.major = element_line(color = "grey20", size = .2),
panel.grid.minor = element_blank(),
plot.title = element_text(
face = "bold", size = 40,
color = "white", hjust = .5, vjust = -7
),
plot.caption = element_text(
size = 14, color = "white",
hjust = .5, vjust = 16
),
plot.subtitle = element_text(
size = 60, color = "#c43c4e",
hjust = .5, vjust = -2
),
plot.margin = unit(c(t = -4, r = -4, b = -4, l = -4), "lines"),
plot.background = element_rect(fill = "grey20", color = NA),
panel.background = element_rect(fill = "grey20", color = NA),
legend.background = element_rect(fill = "grey20", color = NA),
panel.border = element_blank()
) +
labs(
x = "",
y = "",
title = "Active fires in Africa adjusted by confidence level",
subtitle = "{as.Date(frame_time)}",
caption = "©2022 Milos Popovic (https://milospopovic.net)
MODIS Collection 61 NRT Hotspot/Active Fire Detections MCD14DL
from NASA FIRMS"
)
return(africa_map)
}Now, this will create a static map of Africa but to turn it into a timelapse map we need to add the transition time and choose our animation effect. In this tutorial, I opted for “sine-in-out”, which creates a flashy appear/disappear effect that resembles the fire, enriched with a fading effect at the beginning and ending of every frame. There are numerous options in gganimate and I encourage you to check the dedicated page.
In the final part, we define the number of frames: the higher the number of frames, the smoother the transition between every frame will be! I chose 80 frames but you can add more. The same applies to the duration and pause. The reason why I used such a long end_pause is that it will freeze the last frame so that the user can easily gauge the final status.
africa_map <- get_africa_map()
africa_map <- africa_map +
transition_time(acq_date) +
enter_fade() +
exit_shrink() +
ease_aes("sine-in-out", interval = .25)
africa_anim <- animate(africa_map,
nframes = 80,
duration = 25,
start_pause = 3,
end_pause = 30,
height = 6,
width = 7.15,
res = 300,
units = "in"
)
anim_save("active_fires_africa2.gif", africa_anim)Aaaaand, here is the final map! 😍
Alright, that’s all ladies and gents! In this tutorial, you learned another trick on using the satellite data and filtering your area of interest. And all this was achieved in R! Feel free to check the full code here, clone the repo and reproduce, reuse and modify the code as you see fit.
Following this tutorial you can choose another continent or a country. NASA has a similar web service page from which you can query active fires for any country in the world! This means that you can tweak my code a bit to create a myriad of active fire maps.
I’d be happy to hear your view on how this map could be improved or extended to other geographic realms. To do so, please follow me on Twitter, Instagram or Facebook! Also, feel free to support my work by buying me a coffee here!