added data cute, time series (#11)

tutorials/r/3-extract-satellite-data-within-boundary.qmd

@@ -1,14 +1,13 @@
 ---
 title: Extract data within a boundary
-author: NOAA CoastWatch
+author: NOAA CoastWatch, NOAA Openscapes
 ---
 
 ::: {.callout-note title="Learning Objectives"}
 
-1. How to access and download sea surface temperature from  NASA Earthdata
-2. How to trim satellite data to specific bounding coordinates
-3. How to cropped satellite data
-4. How to apply shapefiles as masks to satellite data 
+1. How to access and download sea surface temperature from NASA Earthdata
+2. How to apply shapefiles as masks to satellite data
+3. How to compute monthly average sea surface temperature
 :::
 
 
@@ -44,9 +43,6 @@ only takes a moment to set up.
 
 *Note: See the set-up tab (in left nav bar) for instructions on getting set up on your own computer.*
 
-# Extract data within a boundary
-
-
 
 ## Datasets used
 __GHRSST Level 4 AVHRR_OI Global Blended Sea Surface Temperature Analysis (GDS2) from NCEI__    
@@ -72,7 +68,7 @@ For this exercise we will use the Gulf Stream province (ProvCode: GFST)
 library(terra)
 library(earthdatalogin)
 library(sf)
-
+library(ggplot2)
 ```
 
 
@@ -101,12 +97,11 @@ ycoord <- st_coordinates(GFST)[,2]
 
 ```
 
-## Search data with the dataset unique name and coordinates/dates 
+## Search data from NASA Earthdata with the dataset unique name and coordinates/dates 
 
 ```{r}
 # Connect to NASA Earthdata with no credentials
 edl_netrc()
-
 ```
 
 ```{r}
@@ -133,60 +128,90 @@ length(results)
 
 There are `r length(results)` files.
 
-## Crop and plot one image
+
+## Apply shapefiles as mask to satellite data
+
 
 ```{r}
 
 # Select the first result
 ras <- terra::rast(results[1], vsi = TRUE)
 
-# Set boundaries to crop
-e <- ext(c(min(xcoord), max(xcoord),  min(ycoord), max(ycoord)))
+# Extract SST from the multi-layer raster data
+ras_sst <- rc[["analysed_sst"]]
 
-# Crop the result
-rc <- terra::crop(ras, e)
-```
+# Vectorize shapes
+shp <- vect(shapes)
 
-This has the following layers.
-```{r}
+# Get boundaries for GFST
+GFST <- shp[shp$ProvCode == "GFST",]
+
+# Plot the SST data
+plot(ras_sst)
+
+# Plot GFST boundaries from shapefile
+plot(GFST,col='red')
 
-# Examine variable names
-names(rc)
+# Mask SST with the GFST boundaries
+masked_rc <- mask(ras_sst, GFST)
 
-# Examine the raster data
-rc
+# Visualize the SST in GFST Province
+plot(masked_rc)
 ```
-
 
-## Apply a mask
+## Compute monthly average of SST
 
-Apply mask to a terra raster.
+We will construct a data cube to compute monthly average for sea surface
+temperature data within the boundary.
 
-https://rdrr.io/cran/terra/man/mask.html
+To minimize data loading times, the first 10 results, which correspond to 
+approximately two months of data, will be used for this exercise.
 
 
 ```{r}
 
-# Extract SST from the multi-layer raster data
-rc_sst <- rc[["analysed_sst"]]
+# Select the first 10 SST results
+ras_all <- terra::rast(results[c(1:10)], vsi = TRUE)
 
+# Trim the SST data to the boundaries of GFST
+rc_all <- terra::mask(ras_all, GFST)
 
-# Vectorize shapes
-shp <- vect(shapes)
+# SST data
+rc_sst <- rc_all["analysed_sst", ]
 
-# Get boundaries for GFST
-GFST <- shp[shp$ProvCode == "GFST",]
+# Function to convert times to year-month format
+year_month <- function(x) {
+  format(as.Date(time(x), format="%Y-%m-%d"), "%Y-%m")
+}
 
-# Plot the SST data
-plot(rc_sst)
+# Convert time to Year-month format for aggregation
+ym <- year_month(rc_all)
 
-# Show GFST boundaries from shapefile
-plot(GFST,col='red')
+# Compute raster mean grouped by Year-month
+monthly_mean_rast <- terra::tapp(rc_all, ym, fun = mean)
 
-# Mask SST with the GFST boundarires
-masked_rc <- mask(rc_sst, GFST)
+# Compute mean across raster grouped by Year-month
+monthly_means <- global(monthly_mean_rast, fun = mean)
+```
 
-# Visualize the SST in GFST
-plot(masked_rc)
+## Convert raster into data frame 
+```{r}
+
+# Convert raster into data.frame
+monthly_means_df <- as.data.frame(monthly_means)
+
+# Convert year_month to a column
+monthly_means_df$year_month <- sub("X", "", rownames(monthly_means_df))
 ```
+## Plot monthly mean of sea surface temperature within GFST province
+
+```{r}
+
+# Plot monthly mean
+ggplot(data = monthly_means_df, aes(x = year_month, y = mean, group = 1)) +
+  geom_line() +
+  geom_point() +
+  xlab("Year.Month") + 
+  ylab("Mean SST (F)")
 
+```