Hundreds of thousands of people in east Africa have been displaced and hundreds have died as a result of torrential rains which ended a drought but saturated soils and engorged rivers, resulting in extreme flooding in 2018.
This post will explore these events using the R package smapr, which provides access to global satellite-derived soil moisture data collected by the NASA Soil Moisture Active-Passive (SMAP) mission and abstracts away some of the complexity associated with finding, acquiring, and working with the HDF5 files that contain the observations (shout out to Laura DeCicco and Marco Sciaini for reviewing smapr, and Noam Ross for editing in the rOpenSci onboarding process).
We will focus on Somalia and Kenya, two of the hardest hit countries.
We’ll also lean on another rOpenSci package, rnoaa, to link precipitation to soil moisture.
First, let’s get spatial boundaries for the study area:
library(raster)
library(tidyverse)
library(smapr)
library(rworldmap)
library(rnoaa)
library(plotly)
library(rasterVis)
library(animation)
library(patchwork)
library(sf)
worldmap The striping occurs because of the orbit of the SMAP satellite.
If we were interested in some weekly measure of soil moisture, we could simply average across days to get a more continuous picture of soil moisture:
weekly_sm Just like we did for the soil moisture data, we’ll bundle up all of these steps into a helper function:
average_precip %
map(cpc_prcp) %__%
bind_rows %__%
filter(precip __= 0) %__%
group_by(lon, lat) %__%
summarize(precip = mean(precip, na.rm = TRUE)) %__%
ungroup %__%
mutate(lon = ifelse(lon __ 180, lon - 360, lon)) %__%
make_precip_raster
}
Fetching global soil moisture and precipitation from 2015 to present
Now we have two functions average_smap and average_precip that we can use to get global soil moisture and precipitation data for any date range.
Next, we’ll use these functions to get data at weekly intervals from the beginning of the SMAP data archive in 2015 through the end of August 2018.
start_dates What does the relationship between precipitation and soil moisture look like?
soil_moisture %__%
left_join(precip) %__%
ggplot(aes(pr, sm)) +
geom_point(alpha = .02) +
my_theme +
ylab('Mean soil moisture (m^3 water per m^3 soil)') +
xlab('Mean precipitation (mm)') +
scale_y_log10() +
scale_x_log10()
## Joining, by = c(&quot;x&quot;, &quot;y&quot;, &quot;date&quot;)
Notice the band of zero-precipitation points that sit on the y-axis.
Among nonzero precipitation values, it seems like there is a nonlinear relationship between precipitation and soil moisture.
So, now let’s take a look at the spring and summer of 2018, when the flooding was worst in Somalia and Kenya.
soil_moisture %__%
filter(date __ as.Date('2018-01-01')) %__%
ggplot(aes(x=x, y=y, fill=sm)) + geom_raster() +
scale_fill_viridis_c(direction = -1,
'Soil moisture') +
facet_wrap(~date, nrow = 7) +
theme_minimal() +
theme(axis.text = element_blank()) +
geom_sf(data = study_area_sf, inherit.aes = FALSE, fill = NA, size = .2) +
xlab('') +
ylab('')
Soil moisture
precip %__%
filter(date __ as.Date('2018-01-01')) %__%
ggplot(aes(x=x, y=y, fill=pr)) + geom_raster() +
scale_fill_gradient(low = 'grey95',
high = 'dodgerblue',
'Precipitation') +
facet_wrap(~date, nrow = 7) +
theme_minimal() +
theme(axis.text = element_blank()) +
geom_sf(data = study_area_sf, inherit.aes = FALSE, fill = NA, size = .2) +
xlab('') +
ylab('')
Precipitation
To get a sense for how these values compare to the overall distribution of values in space, we can plot the same data, but rather than using the soil moisture and precipitation values to color the map, we can color the map using the empirical cumulative distribution function (CDF).
This will give us values between 0 and 1, which tell us the fraction of values below a particular value.
So for instance if the empirical CDF at a particular value gives us 0.9, then 90% of the observations were less than that value.
One quick thing to notice is that the distribution of values is different for every pixel:
sm_ecdf %
group_by(x, y) %__%
mutate(ecdf = ecdf(sm)(sm)) %__%
ungroup
sm_ecdf %__%
ggplot(aes(sm, ecdf, group = interaction(x, y))) +
geom_line(alpha = .1) +
my_theme +
xlab('Soil moisture (m^3 water per m^3 soil)') +
ylab('Empirical cumulative distribution function')
pr_ecdf %
group_by(x, y) %__%
mutate(ecdf = ecdf(pr)(pr)) %__%
ungroup
pr_ecdf %__%
ggplot(aes(pr, ecdf, group = interaction(x, y))) +
geom_line(alpha = .1) +
my_theme +
xlab('Precipitation (mm)') +
ylab('Empirical cumulative distribution function')
A few things to notice about these empirical cumulative distribution functions:
The distributions of soil moisture and precipitation are quite different (e.g., the distribution of precipitation has a lot of zero values, and a long tail).
The empirical CDFs provide a mapping from the range of data to the interval from zero to one, so that we can visualize how any particular amount of soil moisture or precipitation compares to the full distribution of values for each pixel. Essentially, when the empirical CDF gives a value close to 0, that’s a low value (a small fraction of values are less than or equal to it). If the empirical CDF gives a value close to one, that’s a high value (a large fraction of values are less than or equal to it).
Let’s see what the empirical CDF values look like on a map:
sm_ecdf %__%
filter(date __ as.Date('2018-01-01')) %__%
ggplot(aes(x=x, y=y, fill=ecdf)) + geom_raster() +
scale_fill_viridis_c('Soil moisture\nempirical CDF',
option = 'B', direction = -1) +
facet_wrap(~date, nrow = 7) +
theme_minimal() +
theme(axis.text = element_blank()) +
geom_sf(data = study_area_sf, inherit.aes = FALSE, fill = NA, size = .2) +
xlab('') +
ylab('')
precip %__%
group_by(x, y) %__%
mutate(ecdf = ecdf(pr)(pr)) %__%
filter(date __ as.Date('2018-01-01')) %__%
ggplot(aes(x=x, y=y, fill=ecdf)) + geom_raster() +
scale_fill_viridis_c('Precipitation\nempirical CDF',
option = 'B',
direction = -1) +
facet_wrap(~date, nrow = 7) +
theme_minimal() +
theme(axis.text = element_blank()) +
geom_sf(data = study_area_sf, inherit.aes = FALSE, fill = NA, size = .2) +
xlab('') +
ylab('')
Across much of Somalia and Kenya, soils were exceptionally dry prior to March 2018 (light colors), and with the prolonged levels of high rain beginning in March, soils were much more saturated than usual through mid to late May (dark colors).
The rain throughout April fell on wet soils, which are less able to absorb water, leading to overland flow and flooding.
We can put this all together in a visualization that shows the evolution of soil moisture and precipitation:
Here the soil moisture and precipitation values are mapped over time, with line plots below to show the median value (black line), interquartile range (dark grey ribbon), 10% and 90% quantiles (light grey ribbon).
The exceptional nature of the flooding in 2018 shows up as a high and broad peak in soil moisture and precipitation, and the drought that preceded the flooding in 2017 is visible as a period of low rains and dry soil.
Want to contribute?
If you’re psyched on global soil moisture and want to contribute to smapr, there is room to develop support for more level 2 products.
If you’re not familiar with the different levels of NASA data products: level 0 is the raw instrument data; level 1 is a bit more processed, e.g., the raw data with ancillary information; level 2 has derived geophysical variables; level 3 has these variables mapped on a uniform grid; and level 4 data consist of model output and derived variables.
Currently smapr supports the more processed level 3 science grade and level 4 enhanced value products primarily.
For a starting point, check out: https://github.com/ropensci/smapr/issues/35
Acknowledgements
The smapr package was developed in Earth Lab with help from Matt Oakley who worked with the Earth Lab Analytics Hub.
The idea for the package emerged from a NOAA Data Partnership event in 2016, where we began working with the NASA SMAP data and realized that we had all the makings for an R package.