---
title: "Rasters Part II: Analyzing and obtaining environmental rasters"
author: "ENS-215"
date: "17-Feb-2025"
output:
  html_document:
    df_print: paged
    theme: journal
    toc: TRUE
    toc_float: TRUE
---

<br/>

Let's load in the packages that we will use today.  Note that while we will be using the `terra` package throughout our analysis today we will not load in the **terra** package.  Instead when we want to use a function from the `terra` package we will use the double-colon notation, for example, `terra::plot()` allows us to call the `plot()` function from the `terra` package.  The reason we do this is because a number of the functions in the `terra` package has names that are identical to functions in the `tidyverse`packages. This can create conflicts/issues with your code and we can avoid these issues by calling `terra` functions using the double-colon notation.  

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(tmap)
library(sf)
```

<br/>

## Get elevation rasters in R

Land surface elevation (topographic) data is frequently useful/required for environmental and geoscience projects.  Obtaining research-quality topographic data can be somewhat time consuming and cumbersome, however there are packages that allow you to seamlessly download this data directly into R.  Today you'll be introduced to the `elevatr` package, which allows you to seamlessly download elevation rasters for anywhere in the world.  This package pulls data from several, research-quality datasets.  More information on the package and the underlying datasets is [available here](https://cran.r-project.org/web/packages/elevatr/vignettes/introduction_to_elevatr.html)  


First let's load in the `elevatr` package.  Let's also load in the `tigris` package, which we've used in the past.  The `tigris` package will be handy for getting borders (e.g., US state, US county), that we will use to define the area of interest for which we would like to download elevation rasters.  

Note: If you haven't already done so, then you will need to first install both the `elevatr` and the `tigris` packages.  

```{r message = F}
library(elevatr)
library(tigris)
```

<br/> 

Let's get the county boundaries for NY state and then use `filter()` to get just Schenectady county.  
```{r message = F, results= 'hide'}
NY_counties <- counties(state = "NY",
                        cb = TRUE) 
```

```{r}
county_schenectady <- NY_counties %>% 
  filter(NAME == "Schenectady")
```


<br/> 

Now that we have the area of interest, let's obtain the elevation raster for that area.  We will use the `get_elev_raster()` function from the `elevatr` package.  The first argument in the function is a sf object that represents our area of interest (i.e., Schenectady county).  The second argument `z` species the resolution of the raster data that we would like to download.  A value of 1 specifies the lowest resolution and a value of 14 is the highest resolution available.  Note that downloading a very high resolution raster for a large area will take much longer (and be a much larger file size).  The `clip` argument specifies if your want to clip the raster to the spatial object.  By specifying `clip = "locations"` the elevation raster will be clipped to the boundaries of Schenectady county.  FYI, when you run the code you will see a download progress bar in your console.  

```{r}
raster_elev <- get_elev_raster(county_schenectady,
                               z = 12,
                               clip = "locations")
```


<br/>

Let's quickly plot the data to confirm that we've downloaded what we expected.  
```{r}
terra::plot(raster_elev)
```



Just to highlight how the raster resolution makes a difference, below I am showing the raster for Schenectady County where `z = 6` (low resolution).  You can see it has much less detail than the raster above where `z = 12` (high resolution).  

**Low resolution** 
```{r message = F, echo= F}
raster_elev_05 <- get_elev_raster(county_schenectady,
                               z = 6,
                               clip = "locations")
```
```{r echo = F} 
terra::plot(raster_elev_05) 
```




<br/> 

Pretty cool - we've quickly and seamlessly obtained raster elevation data for our area of interest.  You can try to modify the code above to test out the `get_elev_raster()` function some more and to try it on different areas.  


<br/>

The `get_elev_raster()` can be used for any location around the world.  To demonstrate this let's obtain raster elevation data for Lesotho (country in southern Africa).  

To obtain high resolution border data we will use the package `rnaturalearth`.  This package let's you easily download a range of spatial data into R.  More info on the package can be found here.  FYI, you will need to install `rnaturalearth` if you haven't yet installed it.  

```{r}
library(rnaturalearth)
```


You should also install the `rnaturalearthhires` package.  

> To do this you will need to first run the following line of code **in your console**.

> `install.packages("remotes")`  


> Then you should run the following line of code **in your console**.  

> `remotes::install_github("ropensci/rnaturalearthhires")`  



Now let's download the country border for Lesotho  
```{r}
borders_hires <- ne_countries(country = "Lesotho",
                              scale = "large", 
                              returnclass = "sf")
```


<br/> 

Following the same approach as we did with Schenectady county, we will now get elevation data for the country of Lesotho  
```{r}
raster_elev_lesotho <- get_elev_raster(borders_hires,
                                       z = 9, 
                                       clip = "locations")
```

<br/>

Let's use `tmap` to make a quick map of the topography of Lesotho  
```{r}
raster_elev_lesotho %>% 
  tm_shape() +
  tm_raster(style = "cont", palette = terrain.colors(n = 10)) +
  
  tm_shape(borders_hires) +
  tm_borders()
```


<br/>  

## Get climate rasters in R

In addition to elevation data, we frequently need gridded (raster) climatic/meteorological data for environmental projects.  Once again, R comes to the rescue and allows you to seamlessly query and download a range of research-quality gridded climatic/meteorological datasets.  One package that is particularly useful in this area is `climateR`.  The `climateR` package gives you access to data from a host of datasets including gridMET (US meteorological data) and TerraClim (global meteorological data).  Each of these datasets has their particular set of variables (e.g., precipitation, temperature, barometric pressure,...) and their particular characteristics (e.g., spatial and temporal resolution, underlying data used to generated the gridded data product).  When conducting research you would want to be aware of these technical details and you would choose the best available data for your application.  To learn more about the `climateR` package and the data it can query you can go to the [following website](https://mikejohnson51.github.io/climateR/), which provides a general overview.  


<br/>


Once you have the remotes package installed (which you should have done earlier in today's lecture), you should then install the `climateR` package by running the following code **in your console**.  

> `remotes::install_github("mikejohnson51/climateR")`


<br/>

After the `climateR` package has been installed, you can load it in using the code below.  
```{r message = F}
library(climateR)
```



```{r echo = F, eval = F}
# Let's run this code to change a setting that will allow us to properly query data when using the `climateR` package.  
sf::sf_use_s2(FALSE)
```

<br/> 

Let's use the `ne_states()` function from the `rnaturalearth` package, which will allow us to download high-resolution borders for the US.  You will likely need to install this package first.  

```{r}
border_NY <- ne_states(country = "United States of America", returnclass = "sf")

border_NY <- border_NY %>% 
  filter(name == "New York")
```


### gridMET data (US meteorological data)
gridMET is a gridded meteorological dataset for the United States.  gridMET is a high-quality dataset that is widely used in climate and environmental research applications.  More details on the data are [available here](https://www.climatologylab.org/gridmet.html).  

Using the `climateR` package you can quickly and efficiently load gridMET data right into R.  Let's first see what variables are available for download.  

```{r}
catalog %>% 
  filter(id == "gridmet") %>% 
  select(varname, units, description)
```

<br/>

Let's go ahead and download the gridded data for the daily minimum air temperature for the state of NY.  The gridMET data is available at daily temporal resolution.  

<br/>

To download the data we will use the `climateR` package.  Let's just download one day of data (FYI, you can easily download multiple days if you want to).  Since Feb 4th, 2023 was an extremely cold day in NY, let's get data for that day.  

to get gridMET data we use the aptly named `getGridMET()` function.  The first argument (`AOI = border_NY`) specifies a spatial polygon (e.g., border) that defines the area over which to get the data.  The second argument specifies the gridMET variable(s) to download.  The third argument specifies the date for which to download data for.  FYI, to download data over more than one day, you would add another argument that specifies the `endDate` to use.  

```{r}
p <- getGridMET(AOI = border_NY,
             varname = c("tmmn"), 
             startDate = "2023-02-04"
             )
```

```{r}
p <- p$daily_minimum_temperature
```

<br/>

Since the gridMET data is obtained as tiles (large rectangular area that covers your area of interest), we often want to crop and mask the data to keep just the pixels that fall within our specified area.  Let's go ahead and do this.  



```{r}
p <- terra::crop(p, border_NY)

p <- terra::mask(p, border_NY)
```

```{r}
p
```
You can see above that the units of temperature are degrees Kelvin.  To convert to degrees C we will simply need to subtract 273.15 from the raster.  

```{r}
p <- p - 273.15
```


<br/>

Now that we have our data and have cropped/masked it to our area of interest, let's make a map of the data.  
```{r}
p %>% 
  tm_shape() + 
  tm_raster(palette = "Blues", n = 10) +
  
  tm_shape(border_NY) +
  tm_borders()
```

<br/>



### TerraClim (Global meteorological data)

gridMET provides excellent daily meteorological data for the US, however if we want gridded met data for outside of the US we will need to rely on another data product.  TerraClim provides two sets of gridded climate data namely, climate normals and climate monthlies.  

Climate normals are gridded monthly average values for a given time period (e.g., 1981-2010).  Thus the January normal of precipitation would be a gridded raster of average precipitation where the average was taken across all of the Januaries between 1981-2010.  

The climate monthlies are gridded monthly data for a given year-month.  Thus you could obtain the gridded precipitation data for January 2020 (or February 2020, or March 1995,...).  


Let's go ahead and get the climate normals for precipitation for the country of Lesotho.  

First, we will use the `ne_countries()` function from the `rnaturalearth` package to get the border of the country of Lesotho.  

```{r}
borders_hires <- rnaturalearth::ne_countries(country = "Lesotho", 
                                             scale = "large", 
                                             returnclass = "sf")
```

<br/>

The variables available from the TerraClim normals are the following:
```{r}
catalog %>% 
  filter(id == "terraclim_normals") %>% 
  select(varname, units, description)
```


<br/> 

Now we will use the `getTerraClimNormals()` from `climateR` to get the TerraClim precipitation normals.  The `AOI` specifies  the spatial object that describes the area of interest.  The `varname` argument specifies the variable dowload.  The `scenario` argument specifies the 30 year normal period to use.  The `month` arugment specifies which months to obtain.  Below we specified `month = 1:12` so we will download all twelve months of average precipitation (averaged over the time period of 1981-2010).

```{r}
climate_raster <- getTerraClimNormals(
  AOI = borders_hires,
  varname = "ppt",
  scenario = "19812010",
  month = 1:12,
  verbose = T,
  dryrun = FALSE
)
```



<br/>

Since `climateR` returns the data as a list of rasters, we will run the following code to obtain the raster from the list object.  This just makes the data easier to deal with in the following steps, as we now have `climate_raster` as a raster object.  
```{r}
climate_raster <- climate_raster$ppt
```


<br/>

Now let's crop and mask the raster.
```{r}
climate_raster <- terra::crop(climate_raster, borders_hires)

climate_raster <- terra::mask(climate_raster, borders_hires)
```



And then make a map.  
```{r}
climate_raster %>% 
  tm_shape() +
  tm_raster(palette = "RdBu", n = 10) +
  
  tm_shape(borders_hires) +
  tm_borders()
  
```

Pretty cool!  We have a faceted map showing the monthly average precipitation (for the period of 1981-2010) for the country of Lesotho.  You could easily modify the above code to obtain precipitation (or other meteorological data) for a different area of interest.  

<br/>

Recalling what we learned in a past class lesson, we know that we can perform mathematical operations on rasters.  It would be useful to sum up the monthly precipitation at each pixel of the raster to obtain the annual average precipitation for each pixel.  This would then give us a raster of the annual average precipitation for Lesotho.  


Let's create this new raster by taking the sum across the twelve monthly rasters.  

```{r}
lesotho_tot_prcp <- sum(climate_raster)
```


<br/>

Now let's make a map to confirm that we did what we had expected to do.  Note that within the `lesotho_tot_prcp` object, the default layer name is `sum`, and thus we need to refer to that layer in the code below.  

```{r}
lesotho_tot_prcp$sum %>% 
  tm_shape() +
  tm_raster(palette = "RdBu", n = 10) + 
  
  tm_shape(borders_hires) +
  tm_borders()
  
```


<br/> 

Now let's use the `climateR` package to get  climate monthlies for a given year-month.  Let's obtain the gridded total precipitation for January 2021 for the country of Lesotho.

```{r}
climate_raster_monthly <- getTerraClim(AOI = borders_hires, 
                                       varname = "ppt", 
                                       startDate = "2021-01-01")
```

```{r}
climate_raster_monthly <- climate_raster_monthly$ppt
```


<br/>

Now let's crop and mask the raster.
```{r}
climate_raster_monthly <- terra::crop(climate_raster_monthly, borders_hires)

climate_raster_monthly <- terra::mask(climate_raster_monthly, borders_hires)
```

<br/>


Now let's make a map of the precipitation data for Jan 2021.  
```{r}
climate_raster_monthly %>% 
  tm_shape() +
  tm_raster(palette = "RdBu", n = 10) + 
  
  tm_shape(borders_hires) +
  tm_borders()
```

## Exercises

Explore the `elevatr` and `climateR` packages.  Try to download additional variables for different areas of interest.  Make some new/modified meteorological and/or topographic maps.  Come up with some interesting questions that you can use these packages to help you answer.  


## Additional rasters  

We will see some additional topics in spatial and raster analysis in the coming lectures, though this is a large topic and we will just scratch the surface in this class.  Let me know if you want to try something else but aren't sure how to move forward.  I can also point you to additional data.  A great place to quickly find some interesting rasters is at [NASA Earth Observation website](https://neo.gsfc.nasa.gov/dataset_index.php).  You can quickly view rasters and download data for preliminary analysis at this site.  The site also provides information and links to the research quality versions of the data if you are interested (though for most initial analysis the data on the link I provide are sufficient).  

```{r eval = F, echo= F}
loc_df <- data.frame(long = -73.902778, lat = 44.365833)
loc_df <- st_as_sf(x = loc_df, coords = c("long", "lat"))

st_crs(loc_df) <- 4326

loc_df <- st_buffer(loc_df, dist = 7500)

raster_elev_mt <- get_elev_raster(loc_df, z = 12)
```

```{r echo = F, eval = F}
terra::plot(raster_elev_mt)
```

```{r echo = F, eval = F}
r_slope <- terra::terrain(x = raster_elev_mt, "slope", unit = "degrees")

terra::plot(r_slope)
```
```{r eval = F, echo = F}
r_TRI <- terra::terrain(x = raster_elev_mt, "TRI", unit = "degrees")

terra::plot(r_TRI)
```



```{r eval = F, echo = F}
## Basemaps
library(basemaps)
```


```{r eval = F, echo = F}
loc_path <- basemap_geotif(border_NY, 
                           map_service = "esri", 
                           map_type = "world_imagery",
                           map_res = 1.0)

loc_path
```






```{r eval = F, echo = F}
loc_basemap <- raster::stack(loc_path)

loc_basemap

#raster::plotRGB(loc_basemap)
```

```{r eval = F, echo = F}

map_full <- tm_shape(loc_basemap) + 
  tm_rgb() +
  tm_shape(border_NY) +
  tm_borders(col = "red") 

map_full
```

```{r eval = F, echo = F}
loc_path <- basemap_geotif(county_schenectady, 
                           map_service = "osm", 
                           map_type = "streets",
                           map_res = 1.0)

loc_path
```

```{r eval = F, echo = F}
loc_basemap <- raster::stack(loc_path)

loc_basemap
```

```{r eval = F, echo = F}
map_full <- tm_shape(loc_basemap) + 
  tm_rgb() + 
  tm_shape(county_schenectady) + 
  tm_borders(col = "red") 

map_full
```


```{r eval = F, echo = F}
## Terrain analysis: Schenectady example
I want to introduce you to the `terrain()` function from the raster package.  This function computes terrain characteristics on an elevation raster.  In particular, you can use the `terrain()` function to compute the slope and aspect (i.e., direction a sloped surface faces) at each pixel.  The `terrain` function also allows for the calculation of a number of other terrain characteristics (see the help file for additional info).  


First let's load in a raster with elevation data for Schenectady.  
```



```{r eval = F, echo = F}
elevation <- raster::raster("https://github.com/stahlm/stahlm.github.io/raw/master/ENS_215/Winter_2022/Lectures/Data/elevation_schdy.tif")
```

<br/>

```{r eval = F, echo = F}
Slope and aspect have a wide range of hydrologic, ecologic, and engineering applications, so let's go ahead and calculate these values for the Schenectady elevation data.  


To calculate the slope, we pass the `elevation` raster to the `terrain()` function.  You'll see that to compute the slope we set the `opt = "slope"`.  Also note that we've set the units to degrees, so that the values returned for the slope are in degrees.  
```


```{r eval = F, echo = F}

r_schdy_slope <- raster::terrain(elevation, opt = "slope", unit = "degrees")
```


```{r eval = F, echo = F}
Now let's compute the aspect.  We do this by setting `opt = "aspect"`  
```

```{r eval = F, echo = F}
r_schdy_aspect <- raster::terrain(elevation, opt = "aspect", unit = "degrees")
```


```{r eval = F, echo = F}
Let's  plot the slope and aspect rasters to see what they look like
```

```{r eval = F, echo = F}
tm_shape(r_schdy_slope) +
  tm_raster(style = "quantile", n = 10) +
  tm_layout(legend.outside = T)
```

<br/>


```{r eval = F, echo = F}
tm_shape(r_schdy_aspect) +
  tm_raster(n = 10) +
  tm_layout(legend.outside = T)
```


<br/>


```{r eval = F, echo = F}
Note that an aspect of 0 degrees indicates that the surface faces north, an aspect of 90 degrees is east, 180 is west,...   

Since the aspect is the direction that a sloping face is pointed, the aspect is less meaningful in relatively flat areas.  Let's create a new raster that has values equal to the aspect when the slope is above a threshold of 1 degree and has a value of `NA` when the slope is less than 1 degree.  

To do this we will first create a raster that has a binary flag (1 if slope >= 1 degree, 0 if slope < 1 degree).  Recall how we replaced raster values in the world elevation raster at the beginning of the lecture.  We'll apply that same approach to create a raster with this binary flag. 
```
 

```{r eval = F, echo = F}
r_schdy_slope_flag <- r_schdy_slope

r_schdy_slope_flag[r_schdy_slope_flag < 1] <- NA
r_schdy_slope_flag[r_schdy_slope_flag >= 1] <- 1


```


```{r eval = F, echo = F} 
Let's take a look at the binary flag raster that we've created
```


```{r eval = F, echo = F}
tm_shape(r_schdy_slope_flag) +
  tm_raster(colorNA = "grey") +
  tm_layout(legend.outside = T)
```


```{r eval = F, echo = F}
You can see that the areas with a slope < 1 degree (flat areas) have NA values (show up grey) and areas with slopes greater than 1 degree all have a value of 1 (show up as tan).  

Now let's create a raster that has the aspect for all pixels with slopes > 1 degree and has values of NA for all other pixels.  This is easy to do since we can simply multiply our `r_schdy_aspect` raster by the `r_schdy_slope_flag` raster.  
```


```{r eval = F, echo = F}
r_schdy_aspect_steep <- r_schdy_slope_flag*r_schdy_aspect
```


<br/> 


```{r eval = F, echo = F}
tm_shape(r_schdy_aspect_steep) +
  tm_raster(n = 12, colorNA = "grey") +
  tm_layout(legend.outside = T)
```




<br/>






 
```{r echo=FALSE, eval=FALSE}
#https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD_LSTD_CLIM_M
```


```{r echo=FALSE, eval=FALSE}
#load(url(("https://github.com/stahlm/stahlm.github.io/raw/master/ENS_215/Winter_2020/Lectures/Spatial_data/MODIS_NPP_2016_07.TIFF")))
```


```{r echo=FALSE, eval=FALSE}
#load(url("https://github.com/stahlm/stahlm.github.io/raw/master/ENS_215/Winter_2020/Lectures/testa.RData"))
```


```{r echo=FALSE, eval=FALSE}
raster_npp_july <- raster::raster("./Spatial_data/MODIS_NPP_2016_07.TIFF")
raster_npp_jan <- raster::raster("./Spatial_data/MODIS_NPP_2016_01.TIFF")
```

```{r echo=FALSE, eval=FALSE}
raster_npp_july
raster_npp_jan
```


```{r echo=FALSE, eval=FALSE}
raster_npp_jan[raster_npp_jan > 50] <- NA
raster_npp_july[raster_npp_july > 50] <-NA
```


```{r echo=FALSE, eval=FALSE}
raster::maxValue(raster_npp_jan)
```


```{r echo=FALSE, eval=FALSE}
raster::plot(raster_npp_jan)
raster::plot(raster_npp_july)
```

```{r echo=FALSE, eval=FALSE}
world_borders <- spData::world
```



```{r echo=FALSE, eval=FALSE}
tm_shape(raster_npp_july) +
  tm_raster() +
  tm_shape(world_borders) +
  tm_borders() +
  tm_layout(legend.position = c("LEFT","BOTTOM"))
```


```{r echo=FALSE, eval=FALSE}
country_npp_july <- raster::extract(raster_npp_july, world_borders, 
                                    fun = mean, na.rm = TRUE, weights = TRUE, sp = TRUE)
```



```{r echo=FALSE, eval=FALSE}
raster_npp_diff_july_jan <- raster_npp_july < raster_npp_jan
```

```{r echo=FALSE, eval=FALSE}
raster::plot(raster_npp_diff_july_jan)
```


```{r echo=FALSE, eval=FALSE}
country_npp_july <- st_as_sf(country_npp_july)
```


```{r echo=FALSE, eval=FALSE}
km2_to_m2 <- 1000*1000

country_npp_july %>% 
  group_by(name_long) %>% 
  summarize(tot_npp = MODIS_NPP_2016_07*area_km2*km2_to_m2) %>% 
  arrange(-tot_npp)
```