Comparing field-collected NPP data to production from the Rangeland Analysis Platform (RAP)
Authors
Greg Maurer
Darren James
Published
July 31, 2026
This tutorial introduces the JRN NPP data and does some comparisons with data from the Rangeland Analysis Platform (RAP). Originally presented in the 2026 Shortcourse workshop.
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.2.1 ✔ readr 2.2.0
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ ggplot2 4.0.3 ✔ tibble 3.3.1
✔ lubridate 1.9.5 ✔ tidyr 1.3.2
✔ purrr 1.2.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(EDIutils)library(sf)
Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
library(terra)
terra 1.9.46
Attaching package: 'terra'
The following object is masked from 'package:tidyr':
extract
library(rapr)library(tidyterra)
Attaching package: 'tidyterra'
The following object is masked from 'package:stats':
filter
# Data access at EDI requires authentication. The easiest way is to create# an API key at https://auth.edirepository.org. Once you have that, place # it in your script. Below, we're getting it from an environmental variable,# but you can replace the "Sys.getenv" part with your actual key.mykey <-Sys.getenv("EDI_API_KEY")# Log in to EDI rogrammatically with an API keyEDIutils::login(key = mykey)
Logged in with EDI-API key.
# Set the dataset IDdatasetID <-"knb-lter-jrn.210011003.106"# Read the list of "entities" in the datasetents <-read_data_entity_names(datasetID)# Read the raw entity dataraw <-read_data_entity(datasetID, ents[1, "entityId"])# Now read the raw data as a CSVanpp_data <- readr::read_csv(file = raw)
Rows: 495 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): zone, site
dbl (2): year, npp_g_m2
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Explore the data a littlehead(anpp_data)
# A tibble: 6 × 4
year zone site npp_g_m2
<dbl> <chr> <chr> <dbl>
1 1990 C CALI 28.6
2 1990 C GRAV 85.7
3 1990 C SAND 103.
4 1990 G BASN 77.2
5 1990 G IBPE 34.1
6 1990 G SUMM 53.7
# Now make a simple plotg <-ggplot(data=anpp_data, aes(x=year, y=npp_g_m2, col=site)) +geom_line()# Show itg
# Now lets subset and plot site==IBPEibpe_anpp <- anpp_data |>filter(site=="IBPE")# Check avialble yearsrange(ibpe_anpp$year) # 1990 to 2022
[1] 1990 2022
# Quick graph of the dataggplot(ibpe_anpp, aes(x = year, y = npp_g_m2)) +geom_line()
# Import shapefile of NPP sites# This can be downloaded from the Jornada Geoportal if you already# have accessNPP_sf <-st_read("../../data/Jornada_prj011_polygons.shp")
Reading layer `Jornada_prj011_polygons' from data source
`/home/runner/work/jeds/jeds/data/Jornada_prj011_polygons.shp'
using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -11895750 ymin: 3827727 xmax: -11878990 ymax: 3851369
Projected CRS: WGS 84 / Pseudo-Mercator
# Look at the columns of the data in the shapefileglimpse(NPP_sf)
Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
# Quick figure of IBPE site polygonggplot(IBPE_sf) +geom_sf()
# Use rapr::get_rap() function to access RAP data# Note that our shapefile is already in the "ESPG:4326" coordinate reference system?get_rap# Currently, RAP biomass and production data are only available as a 30m pixel product# Each pixel is approximately 30m X 3m = 900m^2 in sizest_area(IBPE_sf) /900# The IBPE site is about 5.4 RAP pixels in size
5.429742 [m^2]
# Buffer IBPE area by 30m to ensure complete coverage# Note that the LENGTHUNIT is already in metersIBPE_buffer <- IBPE_sf %>%st_buffer(dist =30)# Figure of IBPE site with 30m bufferggplot() +geom_sf(data = IBPE_buffer) +geom_sf(data = IBPE_sf)
# Use the default terra::plot() function to explore the data# Plots the first 16 rastersplot(IBPE_prod30m)
# Graph the second raster with ggplot and geom_spatraster()# Overlay the IBPE site polygonggplot() +geom_spatraster(data = IBPE_prod30m[[2]]) +geom_sf(data = IBPE_sf, col ="black", fill =NA) +labs(title ="Single Layer from SpatRaster")
# Use the terra::extract() function to calculate mean biomass at the IBPE site for each raster# Using exact = TRUE will caclualte a weigthed mean that uses partial pixels in the site boundary# Using exact = FALSE (default) will only use the pixels whose centroids are within the site boundary (4 pixels in this case)rap_biomass_extract <- terra::extract(IBPE_prod30m, IBPE_sf, fun ="mean", exact =TRUE) # Check the names of the extracted valuesnames(rap_biomass_extract)
# Use tidyr::pivot_longer() to transpose the date from wide to longrap_biomass_df <- rap_biomass_extract %>%pivot_longer(cols =starts_with("vegetation"),names_to ="raster",values_to ="biomass")# Check the text strings of raster names# We need to find a way to extract the year from each of thes text stringsunique(rap_biomass_df$raster)
# We can use the paste() function to construct a vector of all our years# separated by the "bar" symbol which functions as a n "OR" operatorpaste(1990:2022, collapse ="|")
# Use stringr::str_extract() to extract the year from the raster name, then covert to numericrap_biomass_df <- rap_biomass_extract %>%pivot_longer(cols =starts_with("vegetation"),names_to ="raster",values_to ="biomass") %>%mutate(year =str_extract(string = raster, pattern =paste(1990:2022, collapse ="|")) %>%as.numeric())# For each year, sum the biomass for annual forbs and grasses and perennial forbs and grassesrap_biomass_year <- rap_biomass_df %>%group_by(year) %>%summarise(rap_herb_biomass =sum(biomass))# Quick plot of the dataggplot(rap_biomass_year, aes(x = year, y = rap_herb_biomass)) +geom_line()
# We want to join this data (by year) to the filed-based ANPP dataibpe_anpp %>%left_join(rap_biomass_year)
Joining with `by = join_by(year)`
# A tibble: 33 × 5
year zone site npp_g_m2 rap_herb_biomass
<dbl> <chr> <chr> <dbl> <dbl>
1 1990 G IBPE 34.1 540.
2 1991 G IBPE 96.1 827.
3 1992 G IBPE 127. 642.
4 1993 G IBPE 85.7 625.
5 1994 G IBPE 26.6 528.
6 1995 G IBPE 119. 512.
7 1996 G IBPE 97.4 655.
8 1997 G IBPE 137. 671.
9 1998 G IBPE 37.2 487.
10 1999 G IBPE 128 765.
# ℹ 23 more rows
# Convert RAP from lbs/acre to grams/m^2rap_biomass_convert <- rap_biomass_year %>%mutate(rap_g_m2 = rap_herb_biomass *453.59237/4046.85642)# Now merge and transpose from wide to longibpe_biomass_merge <- ibpe_anpp %>%left_join(rap_biomass_convert %>% dplyr::select(year, rap_g_m2)) %>%pivot_longer(cols =c(npp_g_m2, rap_g_m2), names_to ="production_type", values_to ="production_g_m2")
Joining with `by = join_by(year)`
# Quick graph of the dataggplot(ibpe_biomass_merge, aes(x = year, y = production_g_m2, col = production_type)) +geom_line()
# Make the graph prettieribpe_biomass_figure <- ibpe_biomass_merge %>%mutate(source =case_when(production_type =="npp_g_m2"~"Field-based total ANPP (all vegetation except Yucca elata", production_type =="rap_g_m2"~"RAP-based herbaceous biomass (forbs and grasses only)")) # Check the range of the production valuesrange(ibpe_biomass_figure$production_g_m2)
[1] 22.08645 326.60000
ggplot(data = ibpe_biomass_figure, aes(x = year, y = production_g_m2, col = source)) +theme_bw() +geom_line(linewidth =1) +scale_x_continuous(limits =c(1989.5, 2022.5), expand =c(0, 0),breaks =seq(from =1990, to =2022, by =5),minor_breaks =1990:2022) +theme(legend.position ="bottom", legend.direction ="vertical") +xlab("") +ylab(bquote("Production (g*"* m^-2*")")) +ggtitle("Comparison of RAP and field-based production at IBPE site",subtitle ="Years: 1990-2022") +theme(legend.title =element_text(size =14),legend.text =element_text(size =13)) +theme(axis.text =element_text(size =11)) +theme(axis.title =element_text(size =13)) +scale_color_manual(values =c("turquoise", "darksalmon")) +scale_y_continuous(limits =c(0, 340), expand =c(0,0))
