knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=8, fig.height=6 ) NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true") knitr::opts_chunk$set(purl = NOT_CRAN)
library(smapr) library(sp) library(raster)
This vignette outlines a basic use scenario for smapr. We will acquire and process NASA (Soil Moisture Active-Passive) SMAP data, and generate some simple visualizations.
Multiple SMAP data products are provided by the NSIDC, and these products vary in the amount of processing. Currently, smapr primarily supports level 3 and level 4 data products, which represent global daily composite and global three hourly modeled data products, respectively. NSIDC provides documentation for all SMAP data products on their website, and we provide a summary of data products supported by smapr below.
| Dataset id | Description | Resolution | |------------|-----------------------------------------------------|------------| | SPL2SMAP_S | SMAP/Sentinel-1 Radiometer/Radar Soil Moisture | 3 km | | SPL3FTA | Radar Northern Hemisphere Daily Freeze/Thaw State | 3 km | | SPL3SMA | Radar Global Daily Soil Moisture | 3 km | | SPL3SMP | Radiometer Global Soil Moisture | 36 km | | SPL3SMAP | Radar/Radiometer Global Soil Moisture | 9 km | | SPL4SMAU | Surface/Rootzone Soil Moisture Analysis Update | 9 km | | SPL4SMGP | Surface/Rootzone Soil Moisture Geophysical Data | 9 km | | SPL4SMLM | Surface/Rootzone Soil Moisture Land Model Constants | 9 km | | SPL4CMDL | Carbon Net Ecosystem Exchange | 9 km |
This vignette uses the level 4 SPL4SMAU (Surface/Rootzone Soil Moisture Analysis Update) data product.
NASA requires a username and password from their Earthdata portal to access SMAP data. You can get these credentials here: https://earthdata.nasa.gov/
Once you have your credentials, you can use the set_smap_credentials
function to set them for use by the smapr package:
set_smap_credentials("myusername", "mypassword")
This function saves your credentials for later use unless you use the argument
save = FALSE
.
To find out which SMAP data are available, we'll use the find_smap
function,
which takes a data set ID, date(s) to search, and a dataset version.
The SPL4SMAU data product is on version 3 (see https://nsidc.org/data/SPL4SMAU).
available_data <- find_smap(id = 'SPL4SMAU', dates = '2018-06-01', version = 4)
This returns a data frame, where every row is one data file that is available on NASA's servers.
str(available_data)
To download the data, we can use download_smap
. Note that this may take a
while, depending on the number of files being downloaded, and the speed of your
internet connection.
Because we're downloading multiple files, we will use the
verbose = FALSE
argument to avoid printing excessive output to the console.
local_files <- download_smap(available_data, overwrite = FALSE, verbose = FALSE)
Each file corresponds to different times as indicated by the file names:
local_files$name[1:2]
Each file that we downloaded is an HDF5 file with multiple datasets bundled
together.
To list all of the data in a file we can use list_smap
.
By default, if we give list_smap
a data frame of local files, it will
return a list of data frames.
Because all of these data files are of the same data product, using list_smap
on one file (e.g., the first) will tell us what's available in all of the files:
list_smap(local_files[1, ])
To dig deeper, we can use the all
argument to list_smap
:
list_smap(local_files[1, ], all = TRUE)
Looking at this output, we can conclude that the file contains multiple arrays
(notice the dim
column).
These arrays correspond to things like estimated root zone soil moisture
(/Analysis_Data/sm_rootzone_analysis
), estimated surface soil moisture
(/Analysis_Data/sm_surface_analysis
), and estimated surface temperature
(/Analysis_Data/surface_temp_analysis
).
See https://nsidc.org/data/smap/spl4sm/data-fields#sm_surface_analysis for more
detailed information on what these datasets represent and how they were
generated.
The datasets that we are interested in are spatial grids.
The smapr
package can extract these data into raster
objects with the
extract_smap
function, which takes a dataset name as an argument.
These names are paths that can be generated from the output of list_smap
.
For example, if we want to get rootzone soil moisture, we can see a dataset
with name sm_rootzone_analysis
in group /Analysis_Data
, so that the path
to the dataset is /Analysis_Data/sm_rootzone_analysis
:
sm_raster <- extract_smap(local_files, '/Analysis_Data/sm_rootzone_analysis')
This will extract all of the data in the data frame local_files
, generating
a RasterBrick with one layer per file:
sm_raster
We can visualize each layer in our RasterBrick:
plot(sm_raster)
If you want to crop the data to a study region, you can use the raster::crop
function.
Let's illustrate by focusing on the state of Colorado, which is approximately
rectangular, so that we can define an extent
object using latitude and
longitude values that roughly correspond to the state boundaries.
The raster::extent()
function can be used with any Spatial*
or Raster*
object.
co_extent <- extent(c(-109, -102, 37, 41)) co_extent <- as(co_extent, "SpatialPolygons") sp::proj4string(co_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" co_extent
Now that we have a SpatialPolygons object, we can use this to crop the soil moisture raster to our study region. First, we need to ensure that the projections are the same:
proj_co_extent <- spTransform(co_extent, crs(sm_raster))
Then, we could crop the soil moisture data to our polygon, which will make the extents match.
co_soil_moisture <- crop(sm_raster, proj_co_extent) plot(co_soil_moisture)
We might also have a polygon to use as a mask. For example, to mask out our Colorado polygon, we could do the following (operating on the first layer of the soil moisture raster for simplicity):
plot(mask(sm_raster[[1]], proj_co_extent))
Notice that masking does not crop the raster, it simply sets all values outside of the polygon to NA. In some cases, an inverse mask is useful, where values inside the polygon are set to NA:
plot(mask(sm_raster[[1]], proj_co_extent, inverse = TRUE))
We may want to average soil moisture values across layers of a raster brick.
This can be done with the raster::calc()
function:
mean_sm <- calc(sm_raster, fun = mean) plot(mean_sm, main = 'Mean soil moisture')
Our SPL4SMAU data have estimated surface and rootzone soil moisture layers. If we want to compare these values, we can load the surface soil moisture data, compute the mean value over layers as we did for the rootzone soil moisture raster, and generate a scatterplot.
surface_raster <- extract_smap(local_files, name = '/Analysis_Data/sm_surface_analysis') # compute mean mean_surface_sm <- calc(surface_raster, fun = mean) # compare values plot(values(mean_sm), values(mean_surface_sm), col = 'dodgerblue', cex = .1, xlab = 'Rootzone soil moisture', ylab = 'Surface soil moisture', bty = 'n') abline(0, 1, lty = 2)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.