class: left, middle, inverse background-image: url("https://live.staticflickr.com/65535/50362989122_a8ee154fea_k_d.jpg") background-size: cover # .green[Rasters] ### .fancy[Continuously Distributed Data <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:limegreen;overflow:visible;position:relative;"><path d="M149.333 56v80c0 13.255-10.745 24-24 24H24c-13.255 0-24-10.745-24-24V56c0-13.255 10.745-24 24-24h101.333c13.255 0 24 10.745 24 24zm181.334 240v-80c0-13.255-10.745-24-24-24H205.333c-13.255 0-24 10.745-24 24v80c0 13.255 10.745 24 24 24h101.333c13.256 0 24.001-10.745 24.001-24zm32-240v80c0 13.255 10.745 24 24 24H488c13.255 0 24-10.745 24-24V56c0-13.255-10.745-24-24-24H386.667c-13.255 0-24 10.745-24 24zm-32 80V56c0-13.255-10.745-24-24-24H205.333c-13.255 0-24 10.745-24 24v80c0 13.255 10.745 24 24 24h101.333c13.256 0 24.001-10.745 24.001-24zm-205.334 56H24c-13.255 0-24 10.745-24 24v80c0 13.255 10.745 24 24 24h101.333c13.255 0 24-10.745 24-24v-80c0-13.255-10.745-24-24-24zM0 376v80c0 13.255 10.745 24 24 24h101.333c13.255 0 24-10.745 24-24v-80c0-13.255-10.745-24-24-24H24c-13.255 0-24 10.745-24 24zm386.667-56H488c13.255 0 24-10.745 24-24v-80c0-13.255-10.745-24-24-24H386.667c-13.255 0-24 10.745-24 24v80c0 13.255 10.745 24 24 24zm0 160H488c13.255 0 24-10.745 24-24v-80c0-13.255-10.745-24-24-24H386.667c-13.255 0-24 10.745-24 24v80c0 13.255 10.745 24 24 24zM181.333 376v80c0 13.255 10.745 24 24 24h101.333c13.255 0 24-10.745 24-24v-80c0-13.255-10.745-24-24-24H205.333c-13.255 0-24 10.745-24 24z"/></svg>] --- # Rasters .center[ .red[Rasters represent data distributed continuously across a spatial extent] ] -- ### Examples: - Elevation (continuous) - Habitat Type (discrete) - Precipitation (continuous) - Imperviable Surfaces (discrete) --- # What is the Structure of a Raster? A Raster is simply a `matrix` of values with some additional decorations on it that allow it to have a spatial context. ```r values <- rpois( n = 36, lambda=12) values ``` ``` ## [1] 11 13 11 9 10 4 8 19 12 16 14 7 15 10 12 10 8 16 8 15 6 16 11 16 11 ## [26] 16 15 12 9 17 13 11 13 10 9 8 ``` -- ```r x <- matrix( values, nrow=6) x ``` ``` ## [,1] [,2] [,3] [,4] [,5] [,6] ## [1,] 11 8 15 8 11 13 ## [2,] 13 19 10 15 16 11 ## [3,] 11 12 12 6 15 13 ## [4,] 9 16 10 16 12 10 ## [5,] 10 14 8 11 9 9 ## [6,] 4 7 16 16 17 8 ``` --- # Spatial Designations For each value in the `matrix`, when it is turned into a `raster` object: - The `cell` (pixel) has a defined spatial extent (width, height, & origin). - All the physical space represented by that cell has *exactly the same value* - The *courseness* of the raster is question dependent: - 3x5 matrix for Continental US may not adequately capture elevation trends. - 1m² matrix for elevation may be *a bit* too big. --- # Matrix `\(\to\)` Raster ```r library( raster ) r <- raster( x ) r ``` ``` ## class : RasterLayer ## dimensions : 6, 6, 36 (nrow, ncol, ncell) ## resolution : 0.1666667, 0.1666667 (x, y) ## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax) ## crs : NA ## source : memory ## names : layer ## values : 4, 19 (min, max) ``` -- Notice that when I plot it out, it does not show the data, but a summary of the data along with some key data about the contents, including: - A class definition - The dimensions of the underlying data matrix, - The resolution (e.g., the spatial extent of the sides of each pixel). Since we have no CRS here, it is equal to `\(nrows(x)^{-1}\)` and `\(ncols(x)^{-1}\)`. - The extent (the bounding box) and again since we do not have a CRS defined it just goes from `\(0\)` to `\(1\)`. - The `crs` (missing) - The source can be either `memory` if the raster is not that big or `out of memory` if it is just referencing. --- # Loading A Raster By far, you will most commonly working with pre-existing raster data. - Several file formats including GeoTIFF, BIL, & ASC. - All can be loaded from filesystem or internet with address. ```r url <- "https://github.com/dyerlab/ENVS-Lectures/raw/master/data/alt_22.tif" r <- raster( url ) r ``` ``` ## class : RasterLayer ## dimensions : 3600, 3600, 12960000 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -120, -90, 0, 30 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs ## source : https://github.com/dyerlab/ENVS-Lectures/raw/master/data/alt_22.tif ## names : alt_22 ## values : -202, 5469 (min, max) ``` Notice that this raster has a defined CRS and as such it is projected and the extent relates to the units of the datum (e.g., from -120 to -90 degrees longitude and 0 to 30 degrees latitude). --- class: sectionTitle, inverse # .blue[Visualizing Rasters] ## .fancy[What good are they if we cannot look at them?] --- # Built-in Plotting Just like all things in `R`, `raster` objects can be visualized using built-in functions as well as functions from external libraries. .pull-left[ <p>&nbsp;</p> <p>&nbsp;</p> <p>&nbsp;</p> ```r plot( r, xlab="Longitude", ylab="Latitude" ) ``` ] .pull-right[  ] --- # Raster Sizes This particular raster is quite large (in terms of the number of cells) .pull-left[ ``` ## class : RasterLayer ## dimensions : 3600, 3600, 12960000 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -120, -90, 0, 30 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs ## source : alt_22.tif ## names : alt_22 ## values : -202, 5469 (min, max) ``` ] -- .pull-right[ These data only represent the elevation of the land. Where there is water, the value in the underlying matrix is `NA`. Cell Type | Count ----------|-------: Land | 2,469,350 Water | 10,490,650 ] --- # Cropping One of the first things to do is to crop the data down to represent the size and extent of our study area. If we over 10 million missing data points (the ocean) and most of Mexico in this raster above but we are only working with sites in Baja California (Norte y Sur), we would do well to excise (or crop) the raster to only include the area we are interested in working with. -- Workflow: 1. Define a bounding box (the spatial extent of the region of interest) 2. Expand it a bit so that points are not **on the edges** of the box. 3. Create an `extent` 4. Crop the original matrix to represent the boundaries defined in the `extent` --- # Cropping 1: Bounding Box Let's use the beetle data from the [Spatial Points Lecture](https://dyerlab.github.io/ENVS-Lectures/spatial/spatial_points/slides.html#1) as the data we will be working with. ```r library( sf ) library( tidyverse ) beetle_url <- "https://raw.githubusercontent.com/dyerlab/ENVS-Lectures/master/data/Araptus_Disperal_Bias.csv" read_csv( beetle_url ) %>% st_as_sf( coords=c("Longitude","Latitude"), crs=4326 ) -> beetles beetles %>% st_bbox() ``` ``` ## xmin ymin xmax ymax ## -114.29353 23.28550 -109.32700 29.32541 ``` --- # Cropping 2: Expand the Bounding Box .pull-left[ ### Option 1: Eyeball it! ```r beetles %>% st_bbox() ``` ``` ## xmin ymin xmax ymax ## -114.29353 23.28550 -109.32700 29.32541 ``` Maybe rounding it to: ```r eyeball_bbox <- c(-116, -109, 22, 30) ``` ] -- .pull-right[ ### Option 2: Use Buffer ```r beetles %>% st_union() %>% st_buffer( dist = 0.5 ) %>% st_bbox() ``` ``` ## xmin ymin xmax ymax ## -114.29354 23.28549 -109.32699 29.32542 ``` ] --- # Cropping 3: Define the Extent I'll just use the old `eyeball` test to make the numbers 'round'. ```r baja_extent <- extent( eyeball_bbox ) baja_extent ``` ``` ## class : Extent ## xmin : -116 ## xmax : -109 ## ymin : 22 ## ymax : 30 ``` --- # Cropping 4: Cropping To crop the raster, we use the `crop()` function and it makes a new raster (and I can throw the old big one away). ```r alt <- crop( r, baja_extent) rm( r ) # this deletes r from memory alt ``` ``` ## class : RasterLayer ## dimensions : 960, 840, 806400 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -116, -109, 22, 30 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs ## source : memory ## names : alt_22 ## values : -202, 2263 (min, max) ``` --- ```r plot( alt, xlab="Longitude", ylab="Latitude" ) plot( beetles, add=TRUE, col="red", pch=16, cex=1.5) ``` <img src="slides_files/figure-html/unnamed-chunk-13-1.png" width="504" style="display: block; margin: auto;" /> .footnote[(1) Notice the `add=TRUE` adds to the previous plot, and (2) Need to run whole chunk to see built-in plot overlays.] --- # Plotting Rasters with `ggplot` As you probably guessed, there is a `geom_raster()` available to us. .redinline[However], we need to conver the data from a `raster` (`matrix`) to a `data.frame` object that `ggplot` can read. .pull-left[ ```r alt %>% rasterToPoints() %>% head() ``` ``` ## x y alt_22 ## [1,] -115.7958 29.99583 55 ## [2,] -115.7875 29.99583 126 ## [3,] -115.7792 29.99583 94 ## [4,] -115.7708 29.99583 99 ## [5,] -115.7625 29.99583 106 ## [6,] -115.7542 29.99583 120 ``` ] -- .pull-right[ ```r alt %>% rasterToPoints() %>% class() ``` ``` ## [1] "matrix" "array" ``` ] --- # Converting A `raster` to a `data.frame` A little coercion to move `matrix` into `as.data.frame()` is necessary. I also use the `transmute()` function which does in-place renaming (rather than `mutate( X=y ) %>% select( -y )`) ```r alt %>% rasterToPoints() %>% as.data.frame() %>% transmute(Longitude=x, Latitude=y, Elevation=alt_22) -> alt.df head( alt.df ) ``` ``` ## Longitude Latitude Elevation ## 1 -115.7958 29.99583 55 ## 2 -115.7875 29.99583 126 ## 3 -115.7792 29.99583 94 ## 4 -115.7708 29.99583 99 ## 5 -115.7625 29.99583 106 ## 6 -115.7542 29.99583 120 ``` --- # `geom_raster()` .pull-left[ <p>&nbsp;</p> <p>&nbsp;</p> ```r alt.df %>% ggplot() + geom_raster( aes( x = Longitude, y = Latitude, fill = Elevation) ) + coord_equal() + theme_minimal() -> baja_elevation baja_elevation ``` ] .pull-right[ <img src="slides_files/figure-html/unnamed-chunk-18-1.png" width="504" style="display: block; margin: auto;" /> ] --- # Playing with Colors .pull-left[ <img src="slides_files/figure-html/unnamed-chunk-19-1.png" width="504" style="display: block; margin: auto;" /> ] .pull-right[ ### Using Color Gradients There is a built-in `terrain.colors()` function that estimates a set of colors that look somewhat topologically orientated. ```r baja_elevation + scale_fill_gradientn( colors=terrain.colors(100)) ``` ] --- # Custome Color Gradients .pull-left[ Set up a custom gradient with a `low`, `mid`, and `high` color and define the value of the elevation that represents the middle of the range. ```r baja_elevation + scale_fill_gradient2( low = "darkolivegreen", mid = "yellow", high = "brown", midpoint = 1000 ) -> baja_map baja_map ``` ] .pull-right[ <img src="slides_files/figure-html/unnamed-chunk-22-1.png" width="504" style="display: block; margin: auto;" /> ] --- # Overlay Data .pull-left[ <img src="slides_files/figure-html/unnamed-chunk-23-1.png" width="504" style="display: block; margin: auto;" /> ] .pull-right[ Map the `sf` object *over* the background raster and pull it all together. ```r baja_map + geom_sf( aes(size = MFRatio ), data = beetles, color = "dodgerblue2", alpha = 0.75) ``` ] --- class: sectionTitle, inverse # .orange[Raster Manipulations] ## .fancy[Getting Data From the Map] --- # Map Interactivity You can work with raster data interactively (I just cannot do it here on this presentation because it has to be in real time). .pull-left[ ```r plot( alt ) click( alt, xy=TRUE, value=TRUE, n=3 ) -> points ``` ] .pull.right[  ] --- # Points from `click()` Here are what the points look like. ```r points ``` ``` ## x y value ## 1 -113.6292 28.45417 870 ## 2 -112.4792 26.85417 1185 ## 3 -111.2458 24.83750 135 ## 4 -109.9958 23.48750 1145 ``` --- # Reprojecting Rasters Just like points, we can reproject the entire raster using the `projectRaster` function. Here I am going to project the raster into UTM Zone 12N, a common projection for this part of [Mexico from epsg.io](https://epsg.io/6367). Unfortunately, the `raster` library does not use epsg codes so we'll have to use the large description of that projection. See the [page](https://epsg.io/6367) for this projection and scroll down to the proj.4 definition. ```r new.proj <- "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs " ``` --- # Visualizing .pull-left[ ```r alt.utm <- projectRaster( alt, crs = new.proj ) alt.utm ``` ``` ## class : RasterLayer ## dimensions : 981, 879, 862299 (nrow, ncol, ncell) ## resolution : 834, 923 (x, y) ## extent : -21583.09, 711502.9, 2428482, 3333945 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs ## source : memory ## names : alt_22 ## values : -202, 2222.994 (min, max) ``` ] .pull-right[ <img src="slides_files/figure-html/unnamed-chunk-30-1.png" width="504" style="display: block; margin: auto;" /> ] --- # Extracting Data From Rasters > What are the parts of Baja California that are within 100m of the elevation of site named *San Francisquito* (`sfran`)? To answer this, we have the following general outline of operations. 1. Find the coordinates of the site named `sfran` 2. Extract the elevation from the `alt` raster that is within 100m (+/-) of that site. 3. Plot the whole baja data as a background 4. Overlay all the locations within that elevation band. To do this we will use both the `alt` and the `beetles` data objects. --- # Isolating the Point ```r sfran <- beetles$geometry[ beetles$Site == "sfran"] sfran ``` ``` ## Geometry set for 1 feature ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -112.964 ymin: 27.3632 xmax: -112.964 ymax: 27.3632 ## Geodetic CRS: WGS 84 ``` --- # Extracting Data at a Point To extract data from `raster` objects, we need to coerce and specify. ```r raster::extract( alt, as(sfran,"Spatial") ) ``` ``` ## [1] 305 ``` --- .pull-left[ ```r library( ggrepel ) alt.df %>% filter( Elevation >= 205, Elevation <= 405 ) %>% ggplot() + geom_raster( aes( x = Longitude, y = Latitude), fill = "gray80", data = alt.df ) + geom_raster( aes( x = Longitude, y = Latitude, fill = Elevation) ) + scale_fill_gradient2( low = "darkolivegreen", mid = "yellow", high = "brown", midpoint = 305 ) + geom_sf( aes(size=MFRatio), alpha=0.5, color="dodgerblue3", data=beetles) + geom_text_repel( aes( label = Site, geometry = geometry), data = beetles, stat = "sf_coordinates", size = 4, color = "dodgerblue4") + coord_sf() ``` ] .pull-right[ &nbsp;  ] --- class: middle background-image: url("images/contour.png") background-position: right background-size: auto .center[ # Questions?  ] <p>&nbsp;</p> .bottom[ If you have any questions for about the content presented herein, please feel free to [submit them to me](mailto://rjdyer@vcu.edu) and I'll get back to you as soon as possible.]