<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/mapping-the-2018-east-africa-floods-from-space-with-smapr/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ropensci.org/blog/2018/09/25/smapr/"> rOpenSci - open tools for open science</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>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.<br>
This post will explore these events using the R package <code>smapr</code>, which provides access to global satellite-derived soil moisture data collected by the NASA <a rel="nofollow" href="https://smap.jpl.nasa.gov/" target="_blank">Soil Moisture Active-Passive (SMAP) mission</a> and abstracts away some of the complexity associated with finding, acquiring, and working with the <a rel="nofollow" href="https://portal.hdfgroup.org/pages/viewpage.action?pageId=48804287" target="_blank">HDF5</a> files that contain the observations (shout out to <a rel="nofollow" href="https://github.com/ldecicco-USGS" target="_blank">Laura DeCicco</a> and <a rel="nofollow" href="https://github.com/marcosci" target="_blank">Marco Sciaini</a> for reviewing <code>smapr</code>, and <a rel="nofollow" href="https://github.com/noamross" target="_blank">Noam Ross</a> for editing in the rOpenSci onboarding process).<br>
We will focus on Somalia and Kenya, two of the hardest hit countries.<br>
We’ll also lean on another rOpenSci package, <code>rnoaa</code>, to link precipitation to soil moisture.</p>
<p>First, let’s get spatial boundaries for the study area:</p>
<pre><code>library(raster)
library(tidyverse)
library(smapr)
library(rworldmap)
library(rnoaa)
library(plotly)
library(rasterVis)
library(animation)
library(patchwork)
library(sf)

worldmap &lt;- getMap()
study_area &lt;- subset(worldmap, NAME_SORT %in% c('Somalia', 'Kenya'))
study_area_sf &lt;- as(study_area, 'sf')

plot(worldmap)
plot(study_area, add = TRUE, col = 'dodgerblue')
</code></pre>
<p><img alt="" src="https://i2.wp.com/ropensci.org/img/blog-images/2018-09-14-smapr/map-aoi-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>Finding soil moisture data</h2>
<p>Next, we can use the <code>find_smap</code> function to find global soil moisture data for any day.<br>
The SMAP satellite was launched in 2015 with two sensors: an active microwave sensor and a passive radiometer.<br>
The active sensor has since failed, so we will use the radiometer data, specifically version five of the ‘SPL3SMP’ data product (for a full list of data products see <a rel="nofollow" href="https://smap.jpl.nasa.gov/data/" target="_blank">https://smap.jpl.nasa.gov/data/</a>).</p>
<pre><code>find_smap('SPL3SMP', dates = '2018-03-01', version = 5)
</code></pre>
<pre><code>##                               name       date                     dir
## 1 SMAP_L3_SM_P_20180301_R16010_001 2018-03-01 SPL3SMP.005/2018.03.01/
</code></pre>
<p>This returns a data frame with one row per file – we can see here that there is one file available for that date.<br>
If we wanted to search over a range of dates, we could provide a date sequence:</p>
<pre><code>date_seq &lt;- seq(as.Date('2018-03-01'), as.Date('2018-03-06'), by = 1)
files &lt;- find_smap('SPL3SMP', dates = date_seq, version = 5)
files
</code></pre>
<pre><code>##                               name       date                     dir
## 1 SMAP_L3_SM_P_20180301_R16010_001 2018-03-01 SPL3SMP.005/2018.03.01/
## 2 SMAP_L3_SM_P_20180302_R16010_001 2018-03-02 SPL3SMP.005/2018.03.02/
## 3 SMAP_L3_SM_P_20180303_R16010_001 2018-03-03 SPL3SMP.005/2018.03.03/
## 4 SMAP_L3_SM_P_20180304_R16010_001 2018-03-04 SPL3SMP.005/2018.03.04/
## 5 SMAP_L3_SM_P_20180305_R16010_001 2018-03-05 SPL3SMP.005/2018.03.05/
## 6 SMAP_L3_SM_P_20180306_R16010_001 2018-03-06 SPL3SMP.005/2018.03.06/
</code></pre>
<h2>Downloading and extracting soil moisture data</h2>
<p>The <code>download_smap</code> function takes the results from <code>find_smap</code> and downloads files locally:</p>
<pre><code>downloads &lt;- download_smap(files, overwrite = FALSE)
downloads
</code></pre>
<pre><code>##                               name       date                     dir
## 1 SMAP_L3_SM_P_20180301_R16010_001 2018-03-01 SPL3SMP.005/2018.03.01/
## 2 SMAP_L3_SM_P_20180302_R16010_001 2018-03-02 SPL3SMP.005/2018.03.02/
## 3 SMAP_L3_SM_P_20180303_R16010_001 2018-03-03 SPL3SMP.005/2018.03.03/
## 4 SMAP_L3_SM_P_20180304_R16010_001 2018-03-04 SPL3SMP.005/2018.03.04/
## 5 SMAP_L3_SM_P_20180305_R16010_001 2018-03-05 SPL3SMP.005/2018.03.05/
## 6 SMAP_L3_SM_P_20180306_R16010_001 2018-03-06 SPL3SMP.005/2018.03.06/
##                   local_dir
## 1 /home/rstudio/.cache/smap
## 2 /home/rstudio/.cache/smap
## 3 /home/rstudio/.cache/smap
## 4 /home/rstudio/.cache/smap
## 5 /home/rstudio/.cache/smap
## 6 /home/rstudio/.cache/smap
</code></pre>
<p>Next, we can list the contents of the files we downloaded.<br>
These are HDF5 files, so each file contains multiple datasets.<br>
For this data product, a soil moisture data set is named <code>Soil_Moisture_Retrieval_Data_AM/soil_moisture</code>, and we can use the <code>extract_smap</code> function to generate rasters from this dataset (see <code>list_smap</code> for a list of all datasets contained in any file):</p>
<pre><code>sm_raster &lt;- extract_smap(downloads, name = &quot;Soil_Moisture_Retrieval_Data_AM/soil_moisture&quot;)
sm_raster
</code></pre>
<pre><code>## class       : RasterBrick
## dimensions  : 406, 964, 391384, 6  (nrow, ncol, ncell, nlayers)
## resolution  : 36032.22, 36032.22  (x, y)
## extent      : -17367530, 17367530, -7314540, 7314540  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /home/rstudio/.cache/smap/tmp.tif
## names       : SMAP_L3_S//R16010_001, SMAP_L3_S//R16010_001, SMAP_L3_S//R16010_001, SMAP_L3_S//R16010_001, SMAP_L3_S//R16010_001, SMAP_L3_S//R16010_001
## min values  :                  0.02,                  0.02,                  0.02,                  0.02,                  0.02,                  0.02
## max values  :             0.9122642,             0.9390723,             0.9531010,             0.9122642,             0.9531010,             0.9390723
</code></pre>
<p>Let’s plot the RasterBrick to see what the soil moisture data look like:</p>
<pre><code>levelplot(sm_raster)
</code></pre>
<p><img alt="" src="https://i1.wp.com/ropensci.org/img/blog-images/2018-09-14-smapr/plot-smap-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>The striping occurs because of the orbit of the SMAP satellite.<br>
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:</p>
<pre><code>weekly_sm &lt;- mean(sm_raster, na.rm = TRUE)
plot(weekly_sm)
</code></pre>
<p><img alt="" src="https://i2.wp.com/ropensci.org/img/blog-images/2018-09-14-smapr/get-weekly-mean-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<p>We’ll want to do that set of operations a bunch of times for this case study, so we’ll write a little function that takes a date range and returns a raster averaging over dates:</p>
<pre><code>average_smap &lt;- function(date_range) {
  mean_sm &lt;- find_smap('SPL3SMP', dates = date_range, version = 5) %>%
    download_smap(overwrite = FALSE) %>%
    extract_smap(name = &quot;Soil_Moisture_Retrieval_Data_AM/soil_moisture&quot;) %>%
    mean(na.rm = TRUE)
  mean_sm
}
</code></pre>
<p>Notice that spatial coverage is still not 100%.<br>
Notably, places where the satellite missed and frozen regions in the arctic have <code>NA</code> values.</p>
<h2>Getting global precipitation data with rnoaa</h2>
<p>We can better interpret soil moisture if we know about precipitation, and the NOAA Climate Prediction Center (CPC) has global precipitation data that are a cinch to get with the <code>cpc_prcp</code> function.</p>
<pre><code>cpc_prcp('2018-03-01')
</code></pre>
<pre><code>## # A tibble: 259,200 x 3
##      lon   lat precip
##    </code></pre>