solutions 

Rasters

Overview

This notebook will demonstrate how to download Cal-Adapt data as rasters.

Setup

Load caladaptR and the other package we’re going to need. (If you haven’t installed these yet, see this setup script).

library(caladaptr)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
caladaptr (version 0.6.6)
URL: https://ucanr-igis.github.io/caladaptr
Bug reports: https://github.com/ucanr-igis/caladaptr/issues

library(units)

udunits database from C:/Users/Andy/Documents/R/win-library/4.1/units/share/udunits/udunits2.xml

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(lubridate)

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

library(sf)

Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(stars)

Loading required package: abind

library(conflicted)
conflict_prefer("filter", "dplyr", quiet = TRUE)
conflict_prefer("count", "dplyr", quiet = TRUE)
conflict_prefer("select", "dplyr", quiet = TRUE)

Downloading Rasters

The same API Request object can be used to get raster data if you feed it into ca_getrst_stars().

For additional info on downloading and analyzing rasters, see the 3 articles on Downloading Rasters.

Below we get a raster of observed historic temperature data for the Sierra climate region:

sierra_cap <- ca_loc_aoipreset(type = "climregions", idfld = "name", idval = "Sierra") %>% 
  ca_livneh(TRUE) %>% 
  ca_period("year") %>% 
  ca_cvar("pr") %>% 
  ca_years(start = 1970, end = 2010)

sierra_cap

Cal-Adapt API Request
Location(s): 
  AOI Preset: climregions
  name(s): Sierra
Variable(s): pr
Temporal aggregration period(s): year
Livneh data: TRUE
Dates: 1970-01-01 to 2010-12-31
 

sierra_cap %>% ca_preflight(check_for = "getrst")

General issues
 - none found
Issues for downloading rasters
 - none found

plot(sierra_cap, locagrid = TRUE)

To fetch the data as TIFs, use :

tiff_dir <- "./data"

sierra_tiff_fn <- sierra_cap %>% 
  ca_getrst_stars(out_dir = tiff_dir, mask = TRUE, quiet = TRUE, overwrite = FALSE)

Pro Tip:

• to avoid downloading the same TIF multiple times, use the same output directory and set overwrite = FALSE

ca_getrst_stars() returns a vector of TIF files that were downloaded. To work with them, you next have to load them back into R as stars objects (space-time arrays) using ca_stars_read():

sierra_stars_lst <- ca_stars_read(sierra_tiff_fn)
length(sierra_stars_lst)

[1] 1

sierra_stars_lst[[1]]

stars object with 3 dimensions and 1 attribute
attribute(s):
                                 Min.  1st Qu.   Median     Mean  3rd Qu.     Max.   NA's
pr_year_livneh_name-Sierra  0.2009817 2.037323 2.821092 3.106411 3.863086 12.43626 123615
dimension(s):
     from to   offset   delta refsys point values x/y
x       1 58 -121.688  0.0625 WGS 84 FALSE   NULL [x]
y       1 71   40.125 -0.0625 WGS 84 FALSE   NULL [y]
year    1 41     1970       1     NA    NA   NULL    

Plot

To plot a stars objects, you have to decide which layer(s) to plot. In this case, each layer represents a year from 1970 to 2010. Below we plot 4 of the 40 years:

plot(sierra_stars_lst[[1]] %>% slice(index = seq(1,40,length.out =4), along = "year"), 
     axes = TRUE,
     main = attributes(sierra_stars_lst[[1]])$ca_metadata$slug)

Not sure what the units are? You can double-check by viewing the metadata for the slug from the catalog:

ca_catalog_search("pr_year_livneh")

pr_year_livneh
  name: Livneh yearly average precipitation historical
  url: https://api.cal-adapt.org/api/series/pr_year_livneh/
  tres: annual
  begin: 1950-01-01T00:00:00Z
  end: 2013-12-31T00:00:00Z
  units: mm
  num_rast: 64
  id: 382
  xmin: -124.5625
  xmax: -113.375
  ymin: 31.5625
  ymax: 43.75

Subset by Dimension

Subsetting stars objects by a dimension is the counterpart to doing an attribute query in GIS. You can use dplyr filter function as you would a tibble. For example, to get the precipitation for only the 1990s:

sierra_stars_lst[[1]] %>% filter(year >= 1990, year < 2000)

stars object with 3 dimensions and 1 attribute
attribute(s):
                                 Min.  1st Qu.   Median    Mean  3rd Qu.     Max.  NA's
pr_year_livneh_name-Sierra  0.4325525 2.126165 2.905097 3.25017 4.015097 11.33755 30150
dimension(s):

To get the values of a dimension, you can use

(x_range <- sierra_stars_lst[[1]] %>% st_get_dimension_values("x"))

 [1] -121.6562 -121.5938 -121.5312 -121.4688 -121.4062 -121.3438 -121.2812 -121.2188 -121.1562 -121.0938 -121.0312
[12] -120.9688 -120.9062 -120.8438 -120.7812 -120.7188 -120.6562 -120.5938 -120.5312 -120.4688 -120.4062 -120.3438
[23] -120.2812 -120.2188 -120.1562 -120.0938 -120.0312 -119.9688 -119.9062 -119.8438 -119.7812 -119.7188 -119.6562
[34] -119.5938 -119.5312 -119.4688 -119.4062 -119.3438 -119.2812 -119.2188 -119.1562 -119.0938 -119.0312 -118.9688
[45] -118.9062 -118.8438 -118.7812 -118.7188 -118.6562 -118.5938 -118.5312 -118.4688 -118.4062 -118.3438 -118.2812
[56] -118.2188 -118.1562 -118.0938

(yr_range <- sierra_stars_lst[[1]] %>% st_get_dimension_values("year"))

 [1] 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991
[23] 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010

You can also use square bracket notation as you would an array. Just put an extra comma at the beginning of the bracket.

sierra_stars_lst[[1]][ , , , which(yr_range >= 1990 & yr_range < 2000)]

stars object with 3 dimensions and 1 attribute
attribute(s):
                                 Min.  1st Qu.   Median    Mean  3rd Qu.     Max.  NA's
pr_year_livneh_name-Sierra  0.4325525 2.126165 2.905097 3.25017 4.015097 11.33755 30150
dimension(s):

Spatial Selection

To spatially select, you can put a SF object in the first slot of the square bracket notation. Let’s get the precipitation history for Yosemite NP:

(ynp_bnd_sf <- st_read("https://raw.githubusercontent.com/ucanr-igis/caladaptr-res/main/geoms/ynp_bnd.geojson"))

Reading layer `ynp_bnd' from data source 
  `https://raw.githubusercontent.com/ucanr-igis/caladaptr-res/main/geoms/ynp_bnd.geojson' using driver `GeoJSON'
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -119.8864 ymin: 37.4947 xmax: -119.1964 ymax: 38.18515
Geodetic CRS:  WGS 84
Simple feature collection with 1 feature and 11 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -119.8864 ymin: 37.4947 xmax: -119.1964 ymax: 38.18515
Geodetic CRS:  WGS 84
  UNIT_CODE                                                                           GIS_NOTES
1      YOSE Lands - http://landsnet.nps.gov/tractsnet/documents/YOSE/METADATA/yose_metadata.xml
               UNIT_NAME  DATE_EDIT STATE REGION GNIS_ID     UNIT_TYPE CREATED_BY
1 Yosemite National Park 2016-01-27    CA     PW  255923 National Park      Lands
                                                METADATA PARKNAME                       geometry
1 http://nrdata.nps.gov/programs/Lands/YOSE_METADATA.xml Yosemite POLYGON ((-119.8456 37.8327...

plot(ynp_bnd_sf$geometry, axes = TRUE)

To do a spatial selection, put the sf object in the first slot in the square brackets:

ynp_pr_stars <- sierra_stars_lst[[1]][ynp_bnd_sf , , , ]
ynp_pr_stars

stars object with 3 dimensions and 1 attribute
attribute(s):
                                 Min. 1st Qu.  Median    Mean  3rd Qu.     Max. NA's
pr_year_livneh_name-Sierra  0.9827738 2.46228 3.08287 3.32722 3.922521 8.374605 2747
dimension(s):

Plot to verify:

{plot(ynp_pr_stars %>% slice(index = 1, along = "year"), 
     axes = TRUE, main = "Average Daily Precipitation in YNP (mm), 1970", reset = FALSE)
plot(ynp_bnd_sf$geometry, border = "red", lwd = 2, add = TRUE)}

Convert Values / Raster Algebra

Math operators (like * / - +) do ‘pixelwise’ operations on stars objects, just like regular arrays.

To compute the yearly summary in inches, rather than daily average in mm, we just multiply the stars object by a conversion constant. A mm is 0.0393701 inches, so to go from daily averages in mm to annual total in inches:

(ynp_prtotal_stars <- ynp_pr_stars * 0.0393701 * 365)

stars object with 3 dimensions and 1 attribute
attribute(s):
                                Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
pr_year_livneh_name-Sierra  14.12254 35.38317 44.30111 47.81244 56.36696 120.3438 2747
dimension(s):

Plot the first year to verify:

{plot(ynp_prtotal_stars %>% slice(index = 1, along = "year"), 
     axes = TRUE, main = "Total Precipitation in YNP (in), 1970", reset = FALSE)
plot(ynp_bnd_sf$geometry, border = "red", lwd = 2, add = TRUE)}

Aggregate Across Time and Space

To take the average precip per pixel for the entire time range, use st_apply(). The MARGIN argument should be the indices of the dimensions you want to keep:

ynp_prtotal_stars %>% 
  st_apply(MARGIN = 1:2, FUN = mean) %>% 
  plot(axes = TRUE, main = "Average Total Annual Precipitation in YNP (in), 1970-2010")

To take the average total precip of all pixels in the park per year, we simply change the value of MARGIN to 3 (because year is the third dimension). We also have to add na.rm = TRUE to tell the mean function to ignore all the NA values outside the park boundary:

(ynp_pr_annual_avg_stars <- ynp_prtotal_stars %>% 
  st_apply(MARGIN = 3, FUN = mean, na.rm = TRUE))

stars object with 1 dimensions and 1 attribute
attribute(s):
          Min.  1st Qu.  Median     Mean  3rd Qu.     Max.
mean  19.96104 37.36484 41.9291 47.81244 53.77562 90.80345
dimension(s):

stars objects can be converted to data frames:

ynp_pr_annual_avg_df <- ynp_pr_annual_avg_stars %>% as.data.frame
plot(ynp_pr_annual_avg_df, type="b", main = "YNP Total Precip Averaged Across Park (in)")

To plot the distribution of values, we can grab the individual values with [] (just like an array) and make a histogram:

hist(ynp_prtotal_stars, main = "Distribution of Total Annual Precip in YNP")

There is a lot more you can do with rasters, including pixel summaries, combining them into higher dimensional data cubes, spatially mosaicing them, etc. For more info, see the Rasters articles on the website.

Your Turn

Download historic precipitation data for the county where you live or work. [Answer]

## Example: Mendocino County

ca_aoipreset_geom("counties", quiet = TRUE) %>% 
  st_drop_geometry() %>% 
  filter(name == "Mendocino") %>% 
  select(name, state_name, fips)

       name state_name  fips
1 Mendocino California 06045

mendocino_cap <- ca_loc_aoipreset(type = "counties", idfld = "fips", idval = "06045") %>% 
  ca_livneh(TRUE) %>% 
  ca_period("year") %>% 
  ca_cvar("pr") %>% 
  ca_years(start = 1970, end = 2010)

plot(mendocino_cap)

mendocino_fn <- mendocino_cap %>% 
  ca_getrst_stars(out_dir = tempdir(), mask = TRUE, quiet = TRUE, overwrite = FALSE)

mendocino_stars_lst <- ca_stars_read(mendocino_fn) 
  
mendocino_stars_lst[[1]]

stars object with 3 dimensions and 1 attribute
attribute(s):
                                Min.  1st Qu.   Median     Mean  3rd Qu.   Max. NA's
pr_year_livneh_fips-06045  0.9326811 2.792288 3.627976 3.892968 4.744806 11.229 7134
dimension(s):
     from to   offset   delta refsys point values x/y
x       1 20 -124.062  0.0625 WGS 84 FALSE   NULL [x]
y       1 21  40.0625 -0.0625 WGS 84 FALSE   NULL [y]
year    1 41     1970       1     NA    NA   NULL    

plot(mendocino_stars_lst[[1]] %>% slice(index = seq(1,40,length.out =4), along = "year"), 
     axes = TRUE,
     main = attributes(mendocino_stars_lst[[1]])$ca_metadata$slug)

---