Lab07_SpatialMapsPII

Total Points = 100

For this lab activity, you will need: 1. Your final project partner 2. Your research project data!

** ONLY ONE SUBMISSION PER GROUP

All maps should include a: 1. Scale Bar 2. North Arrow 3. Title 4. Legend (if applicable) 5. Appropriate color scale (if applicable) 6. An appropriate CRS

Question 1:

Write a paragraph (4-6 sentences) summarizing your research project in detail. This needs to include your research question. (5)

The purpose of this project is to find correlative relationships between Sea Surface Temperature (SST), Vertical Wind Shear (SHRD), and Rapid Intensification (RI) in the SHIPS model data from 1982–2023. The research question we are answering is “RQ1: Is there a correlation between vertical wind shear and sea surface temperatures in Atlantic Basin tropical cyclones between the years of 1982-2023?” We do this through exploratory data analysis and graphical data analysis with Pearson’s correlation and significance (P-value) testing of our data. Although the question to the thermodynamic relationship between storm intensity, shear, and SST are already proven, this research is evaluating a statistical model that may have different correlative relationships than in theoretical domains. We belive this way, this project adds value to the meteorological field.

Question 2:

Load your data set into R. If possible, make sure your data frame is in an RData file. (5)

AI Used:

“I am having trouble getting the data from output.csv from my .rdata file, what are some ways I can recurrsively search and load the data?”

Output:

” if(file.exists(“your.csv”)) { df <- read.csv(“your.csv”) } else if(file.exists(“yourData.RData”)) { load(“yourData.RData”) df <- get(ls()[1]) } else { df <- read.csv(file.choose()) } ”

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(sf)

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(ggplot2)
library(tmap)
library(USAboundaries)
library(RColorBrewer)
library(stringr)

if(file.exists("output.csv")) {
  storm_raw <- read.csv("output.csv")
} else if(file.exists("finalData.RData")) {
  load("finalData.RData")
  storm_raw <- get(ls()[1]) 
} else {
  storm_raw <- read.csv(file.choose())
}

storm_clean <- storm_raw %>%
  mutate(VMAX = as.numeric(VMAX)) %>%
  group_by(storm) %>%
  arrange(TIME) %>%
  mutate(
    VMAX_future = lead(VMAX, 4),
    DV24 = VMAX_future - VMAX,
    RI_flag = ifelse(!is.na(DV24) & DV24 >= 30, 1, 0)
  ) %>%
  ungroup() %>%
  mutate(
    lat_geo = as.numeric(LAT) / 10,
    lon_geo = -1 * (as.numeric(LON) / 10), 
    SST_Celsius = as.numeric(RSST) / 10,
    Shear_kts = as.numeric(SHRD) / 10,
    TIME = as.numeric(TIME)
  ) %>%
  filter(TIME == 0, !is.na(SST_Celsius), !is.na(Shear_kts), !is.na(lat_geo)) %>%
  mutate(Basin_Region = case_when(
    lon_geo > -100 & lon_geo < -80 & lat_geo > 18 & lat_geo < 31 ~ "Gulf of Mexico",
    lon_geo > -88 & lon_geo < -60 & lat_geo > 9 & lat_geo < 20 ~ "Caribbean Sea",
    TRUE ~ "Atlantic Ocean"
  ))

Question 3:

Spatial visualization is essential in meteorology — from displaying station networks to showing precipitation patterns or model output grids. Therefore, you are tasked with creating a map of your research area that will serve as the foundation for spatial analyses in your final project. For this question, create a base map identifying the region, state, county, city, etc. of your final project. (10)

library(sf)
library(tmap)
library(rnaturalearth) # added from online documentation on cleaner basemap creation and further debugging and integration with AI
library(rnaturalearthdata) # added from online documentation on cleaner basemap creation and further debugging and integration with AI

## 
## Attaching package: 'rnaturalearthdata'

## The following object is masked from 'package:rnaturalearth':
## 
##     countries110

# Get world countries data
world <- ne_countries(scale = 50, returnclass = "sf")

# Create base map using tmap
tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

## ℹ toggle with `tmap::ttm()`
## This message is displayed once per session.

tm_shape(world, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 55), crs = 4326)) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_compass(type = "arrow", size = 3, fontsize = 1, color.dark = "black", position = c("right", "top")) +
  tm_scale_bar(text.size = 0.8, lwd = 3, color.dark = "black", position = c("left", "bottom")) +
  tm_layout(
    title = "Atlantic Hurricane Basin",
    title.size = 1.5,
    title.fontface = "bold",
    title.position = c("center", "top"),
    bg.color = "#e0f3f8",
    inner.margins = c(0.05, 0.05, 0.05, 0.05),
    frame = TRUE
  ) +
  tm_credits("CRS: WGS84. Source: Natural Earth", position = c("center", "bottom"), size = 0.8)

## ! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [tip] Consider a suitable map projection, e.g. by adding `+ tm_crs("auto")`.
## Scale bar set for latitude km and will be different at the top and bottom of the map.

AI used:

“Help me create a bunch of different subplot for the specified basins for the plot I made. I want the plots ot be the zoomed in version of the plot I made (attached). They should each be a new plot with each variable plotted (SST and SHRD)”

Output:

“map_summary_shear <- ggplot() + geom_sf(data = world, fill =”gray90”, color = “gray60”, size = 0.3) + geom_rect(data = gulf_box, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), fill = NA, color = “black”, size = 0.8, linetype = “dashed”) + geom_rect(data = carib_box, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), fill = NA, color = “black”, size = 0.8, linetype = “dashed”) + geom_label(data = basin_summary_data, aes(x = Centroid_Lon, y = Centroid_Lat, label = paste(Basin_Region, “Shear:”, round(Avg_Shear, 1), “kts”)), fill = “white”, alpha = 0.8, fontface = “bold”) + coord_sf(xlim = c(-100, -10), ylim = c(5, 50), crs = 4326, expand = FALSE) + labs(title = “Average Shear by Basin Region”, x = NULL, y = NULL) + theme_bw() + theme( panel.background = element_rect(fill = “#e0f3f8”), panel.grid.major = element_line(color = “gray80”, linetype = “dotted”), axis.text = element_text(size = 12), axis.title = element_text(size = 14), plot.title = element_text(face = “bold”, size = 16, hjust = 0.5), legend.position = “bottom” ) map_summary_shear

create_sst_map(storm_clean, “Gulf: SST”, xlims = c(-98, -80), ylims = c(18, 31)) create_shear_map(storm_clean, “Gulf: Shear”, xlims = c(-98, -80), ylims = c(18, 31)) create_sst_map(storm_clean, “Caribbean: SST”, xlims = c(-88, -60), ylims = c(9, 20)) create_shear_map(storm_clean, “Caribbean: Shear”, xlims = c(-88, -60), ylims = c(9, 20)) create_sst_map(storm_clean, “Atlantic: SST”, xlims = c(-80, -20), ylims = c(10, 45)) create_shear_map(storm_clean, “Atlantic: Shear”, xlims = c(-80, -20), ylims = c(10, 45)) create_sst_map(storm_clean, “Gulf: SST”, xlims = c(-98, -80), ylims = c(18, 31)) create_sst_map(storm_clean, “Caribbean: SST”, xlims = c(-88, -60), ylims = c(9, 20)) create_sst_map(storm_clean, “Atlantic: SST”, xlims = c(-80, -20), ylims = c(10, 45)) create_shear_map(storm_clean, “Gulf: Shear”, xlims = c(-98, -80), ylims = c(18, 31)) create_shear_map(storm_clean, “Caribbean: Shear”, xlims = c(-88, -60), ylims = c(9, 20)) create_shear_map(storm_clean, “Atlantic: Shear”, xlims = c(-80, -20), ylims = c(10, 45)) ”

“── tmap v3 code detected ──────────────────────────────────────────────────────────────────────────────────────── [v3->v4] tm_lines(): use col_alpha instead of alpha. [v3->v4] symbols(): use fill_alpha instead of alpha. Error in tm_scale_continuous(values =”brewer.rd_yl_bu”, direction = -1) : unused argument (direction = -1)

Help me fix this error in my rstudio code chunk”

Output:

“tm_shape(world, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 55), crs = 4326)) + tm_polygons(fill =”gray90”, col = “gray60”, lwd = 0.5) + tm_shape(storm_tracks) + tm_lines(col = “black”, lwd = 0.3, alpha = 0.2) + tm_shape(ri_events_sf) + tm_symbols( col = “SST (°C)”, size = “Shear (Knots)”, col.scale = tm_scale_continuous(values = “brewer.rd_yl_bu”, direction = -1), size.scale = tm_scale(values.range = c(1, 0.3)), shape = 21, alpha = 0.8, col.legend = tm_legend(orientation = “landscape”), size.legend = tm_legend(orientation = “landscape”) ) + tm_shape(world) + tm_polygons(fill = “gray90”, col = “gray60”, lwd = 0.5) + tm_compass(type = “arrow”, size = 3, fontsize = 1, color.dark = “black”, position = c(“right”, “top”)) + tm_scale_bar(text.size = 0.8, lwd = 3, color.dark = “black”, position = c(“left”, “bottom”)) + tm_layout( title = “Locations of Rapid Intensification (RI) Events”, title.size = 1.5, title.fontface = “bold”, title.position = c(“center”, “top”), bg.color = “#e0f3f8”, legend.bg.color = “white”, legend.position = c(“right”, “center”), inner.margins = c(0.05, 0.05, 0.05, 0.05) )”

“I am unable to change the color bar from blue to my blue-red SST color bar. Help me make these changes to my RI plot”

Output ” scale_fill_gradientn(colors = c(“#d7191c”, “#fdae61”, “#abd9e9”, “#2c7bb6”), ”

—this was useless btw did not work for my RI plot.

Question 4:

Using the base map from Question 3, add at least one spatial data layer (points, lines, polygons) of meteorology information. (You can add more than one if appropriate) ** If you are having trouble, consider this an opportunity to explain the importance of your research question. ** (25)

library(sf)
library(tmap)
library(rnaturalearth) # added from online documentation on cleaner basemap creation and further debugging and integration with AI
library(rnaturalearthdata) # added from online documentation on cleaner basemap creation and further debugging and integration with AI
library(ggplot2)
library(dplyr)
library(RColorBrewer)
library(scales)

world <- ne_countries(scale = 50, returnclass = "sf")

ri_events <- storm_clean %>% filter(RI_flag == 1)

ri_events_sf <- st_as_sf(ri_events, coords = c("lon_geo", "lat_geo"), crs = 4326) %>%
  mutate(
    `SST (°C)` = SST_Celsius,
    `Shear (Knots)` = Shear_kts
  )

storm_tracks <- storm_clean %>%
  group_by(storm) %>%
  summarise(
    geometry = st_combine(st_as_sf(cur_data(), coords = c("lon_geo", "lat_geo"), crs = 4326)$geometry),
    .groups = "drop"
  ) %>%
  st_as_sf() %>%
  st_cast("LINESTRING")

## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `geometry = st_combine(...)`.
## ℹ In group 1: `storm = "ALBE"`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.

tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

tm_shape(world, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 55), crs = 4326)) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_shape(storm_tracks) +
  tm_lines(col = "black", lwd = 0.3, alpha = 0.2) +
  tm_shape(ri_events_sf) +
  tm_symbols(
    col = "SST (°C)", 
    size = "Shear (Knots)", 
    col.scale = tm_scale(values = c("#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")),
    size.scale = tm_scale(values.range = c(1, 0.3)),
    shape = 21, 
    alpha = 0.8,
    col.legend = tm_legend(orientation = "landscape"),
    size.legend = tm_legend(orientation = "landscape")
  ) +
  tm_shape(world) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_compass(type = "arrow", size = 3, fontsize = 1, color.dark = "black", position = c("right", "bottom")) +
  tm_scale_bar(text.size = 0.8, lwd = 3, color.dark = "black", position = c("left", "bottom")) +
  tm_layout(
    title = "Locations of Rapid Intensification (RI) Events",
    title.size = 1.5,
    title.fontface = "bold",
    title.position = c("center", "top"),
    bg.color = "#e0f3f8",
    legend.bg.color = "white",
    legend.position = c("right", "center"),
    inner.margins = c(0.05, 0.05, 0.05, 0.05)
  )

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_lines()`: use `col_alpha` instead of `alpha`.[v3->v4] `symbols()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `symbols()`: use `fill_alpha` instead of `alpha`.! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`Scale bar set for latitude km and will be different at the top and bottom of the map.

create_sst_map <- function(data, title_txt, xlims = c(-100, -10), ylims = c(5, 50)) {
  max_sst <- max(data$SST_Celsius, na.rm = TRUE)
  ggplot() +
    stat_summary_2d(data = data, aes(x = lon_geo, y = lat_geo, z = SST_Celsius), 
                    fun = mean, bins = 60, alpha = 0.9) +
    geom_path(data = data, aes(x = lon_geo, y = lat_geo, group = storm), 
              color = "black", size = 0.1, alpha = 0.15) +
    geom_sf(data = world, fill = "gray95", color = "gray70", size = 0.2) +
    scale_fill_gradientn(colors = c("#2c7bb6", "#abd9e9", "#fdae61", "#d7191c"), 
                         name = "SST (°C)", limits = c(24, max_sst), oob = scales::squish) +
    theme_bw() +
    coord_sf(xlim = xlims, ylim = ylims, crs = 4326, expand = FALSE) +
    labs(title = title_txt, subtitle = "Mean Sea Surface Temperature", x = NULL, y = NULL) +
    theme(
      panel.background = element_rect(fill = "#e0f3f8"),
      panel.grid.major = element_line(color = "gray80", linetype = "dotted"),
      axis.text = element_text(size = 12),
      axis.title = element_text(size = 14),
      plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
      plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray40"),
      legend.position = "right"
    )
}

create_shear_map <- function(data, title_txt, xlims = c(-100, -10), ylims = c(5, 50)) {
  ggplot() +
    stat_summary_2d(data = data, aes(x = lon_geo, y = lat_geo, z = Shear_kts), 
                    fun = mean, bins = 60, alpha = 0.9) +
    geom_path(data = data, aes(x = lon_geo, y = lat_geo, group = storm), 
              color = "black", size = 0.1, alpha = 0.15) +
    geom_sf(data = world, fill = "gray95", color = "gray70", size = 0.2) +
    scale_fill_gradientn(colors = c("#d7191c", "#fdae61", "#abd9e9", "#2c7bb6"), 
                         name = "Shear (Knots)") +
    theme_bw() +
    coord_sf(xlim = xlims, ylim = ylims, crs = 4326, expand = FALSE) +
    labs(title = title_txt, subtitle = "Mean Vertical Wind Shear", x = NULL, y = NULL) +
    theme(
      panel.background = element_rect(fill = "#e0f3f8"),
      panel.grid.major = element_line(color = "gray80", linetype = "dotted"),
      axis.text = element_text(size = 12),
      axis.title = element_text(size = 14),
      plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
      plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray40"),
      legend.position = "right"
    )
}

create_sst_map(storm_clean, "Atlantic Basin: SST")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

create_shear_map(storm_clean, "Atlantic Basin: Shear")

basin_summary_data <- storm_clean %>% 
  group_by(Basin_Region) %>% 
  summarise(
    Avg_SST = mean(SST_Celsius, na.rm = TRUE), 
    Avg_Shear = mean(Shear_kts, na.rm = TRUE), 
    Centroid_Lon = mean(lon_geo, na.rm = TRUE), 
    Centroid_Lat = mean(lat_geo, na.rm = TRUE)
  )

gulf_box <- data.frame(xmin = -100, xmax = -80, ymin = 18, ymax = 31)
carib_box <- data.frame(xmin = -88, xmax = -60, ymin = 9, ymax = 20)

map_summary_sst <- ggplot() +
  geom_sf(data = world, fill = "gray90", color = "gray60", size = 0.3) +
  geom_rect(data = gulf_box, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), 
            fill = NA, color = "black", size = 0.8, linetype = "dashed") +
  geom_rect(data = carib_box, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), 
            fill = NA, color = "black", size = 0.8, linetype = "dashed") +
  geom_label(data = basin_summary_data, 
             aes(x = Centroid_Lon, y = Centroid_Lat, 
                 label = paste(Basin_Region, "\nAvg SST:", round(Avg_SST, 1), "°C")), 
             fill = "white", alpha = 0.8, fontface = "bold") +
  coord_sf(xlim = c(-100, -10), ylim = c(5, 50), crs = 4326, expand = FALSE) +
  labs(title = "Average SST by Basin Region", x = NULL, y = NULL) +
  theme_bw() +
  theme(
    panel.background = element_rect(fill = "#e0f3f8"),
    panel.grid.major = element_line(color = "gray80", linetype = "dotted"),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    legend.position = "right"
  )
map_summary_sst

map_summary_shear <- ggplot() +
  geom_sf(data = world, fill = "gray90", color = "gray60", size = 0.3) +
  geom_rect(data = gulf_box, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), 
            fill = NA, color = "black", size = 0.8, linetype = "dashed") +
  geom_rect(data = carib_box, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), 
            fill = NA, color = "black", size = 0.8, linetype = "dashed") +
  geom_label(data = basin_summary_data, 
             aes(x = Centroid_Lon, y = Centroid_Lat, 
                 label = paste(Basin_Region, "\nAvg Shear:", round(Avg_Shear, 1), "kts")), 
             fill = "white", alpha = 0.8, fontface = "bold") +
  coord_sf(xlim = c(-100, -10), ylim = c(5, 50), crs = 4326, expand = FALSE) +
  labs(title = "Average Shear by Basin Region", x = NULL, y = NULL) +
  theme_bw() +
  theme(
    panel.background = element_rect(fill = "#e0f3f8"),
    panel.grid.major = element_line(color = "gray80", linetype = "dotted"),
    axis.text = element_text(size = 12),
    axis.title = element_text(size = 14),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    legend.position = "bottom"
  )
map_summary_shear

create_sst_map(storm_clean, "Gulf: SST", xlims = c(-98, -80), ylims = c(18, 31))

create_shear_map(storm_clean, "Gulf: Shear", xlims = c(-98, -80), ylims = c(18, 31))

create_sst_map(storm_clean, "Caribbean: SST", xlims = c(-88, -60), ylims = c(9, 20))

create_shear_map(storm_clean, "Caribbean: Shear", xlims = c(-88, -60), ylims = c(9, 20))

create_sst_map(storm_clean, "Atlantic: SST", xlims = c(-80, -20), ylims = c(10, 45))

create_shear_map(storm_clean, "Atlantic: Shear", xlims = c(-80, -20), ylims = c(10, 45))

create_sst_map(storm_clean, "Gulf: SST", xlims = c(-98, -80), ylims = c(18, 31))

create_sst_map(storm_clean, "Caribbean: SST", xlims = c(-88, -60), ylims = c(9, 20))

create_sst_map(storm_clean, "Atlantic: SST", xlims = c(-80, -20), ylims = c(10, 45))

create_shear_map(storm_clean, "Gulf: Shear", xlims = c(-98, -80), ylims = c(18, 31))

create_shear_map(storm_clean, "Caribbean: Shear", xlims = c(-88, -60), ylims = c(9, 20))

create_shear_map(storm_clean, "Atlantic: Shear", xlims = c(-80, -20), ylims = c(10, 45))

Question 5:

Time to make it more complicated. Conduct a geocomputation (st_intersects, st_join, st_buffer, st_contains, etc.) on your data. Make a map of the information that you received. (35) Here are some examples: 1. Number of survey respondents in each county/state. 2. Average Height of storm surge in coastal counties 3. Total number of tornadoes by county in the United States 4. Maximum CAPE value found in a single tornado outbreak 5. Total number of weather delays by airport

world = ne_countries(scale = 50, returnclass = "sf")

basin_bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 50), crs = 4326)

grid_poly = st_make_grid(st_as_sfc(basin_bbox), cellsize = 2, square = TRUE)
grid_poly = st_sf(geometry = grid_poly)
grid_poly = grid_poly %>% mutate(grid_id = row_number())

suppressWarnings({
  grid_centroids <- st_centroid(grid_poly)
})

grid_coords <- st_coordinates(grid_centroids)


grid_poly = st_make_grid(st_as_sfc(basin_bbox), cellsize = 2, square = TRUE)
grid_poly = st_sf(geometry = grid_poly)
grid_poly = grid_poly %>% mutate(grid_id = row_number())
grid_coords <- st_coordinates(grid_centroids)
grid_poly = grid_poly %>%
  mutate(
    lon = grid_coords[,1],
    lat = grid_coords[,2],
    Basin_Region = case_when(
      lon > -100 & lon < -80 & lat > 18 & lat < 31 ~ "Gulf of Mexico",
      lon > -88 & lon < -60 & lat > 9 & lat < 20 ~ "Caribbean Sea",
      TRUE ~ "Atlantic Ocean"
    )
  )

basin_lookup <- grid_poly %>%
  st_drop_geometry() %>%
  select(grid_id, Basin_Region)

storm_points_sf <- st_as_sf(storm_clean, coords = c("lon_geo", "lat_geo"), crs = 4326)

joined_data <- st_join(grid_poly, storm_points_sf, join = st_intersects)

ri_probability <- joined_data %>%
  group_by(grid_id) %>%
  summarise(
    Total_Points = n(), 
    RI_Events = sum(RI_flag, na.rm = TRUE), 
    Probability = (RI_Events / Total_Points) * 100,
    .groups = "drop"
  ) %>%
  filter(Total_Points > 10)

tmap_mode("plot")

## ℹ tmap modes "plot" - "view"

tm_shape(ri_probability, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 50), crs = 4326)) +
  tm_polygons(
    fill = "Probability", 
    fill.scale = tm_scale(values = c("white", "#d7191c")),
    fill.legend = tm_legend(title = "RI Probability (%)", orientation = "portrait", position = tm_pos_in("right", "center")), 
    col = "white", 
    lwd = 0.1
  ) +
  tm_shape(world) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_compass(type = "arrow", size = 3, fontsize = 1, color.dark = "black", position = c("right", "top")) +
  tm_scalebar(text.size = 0.8, lwd = 3, color.dark = "black", position = c("left", "bottom")) +
  tm_title(text = "Probability of Rapid Intensification by Region", size = 1.5, fontface = "bold", position = c("center", "top")) +
  tm_layout(
    bg.color = "#e0f3f8", 
    legend.bg.color = "white",
    inner.margins = c(0.1, 0.05, 0.05, 0.05)
  )

## Scale bar set for latitude km and will be different at the top and bottom of the map.

p_value_stats <- joined_data %>%
  group_by(grid_id) %>%
  summarise(
    Count = n(),
    P_Value = if(n() > 5 && sd(SST_Celsius, na.rm=TRUE) > 0 && sd(Shear_kts, na.rm=TRUE) > 0) {
      cor.test(SST_Celsius, Shear_kts, method = "pearson")$p.value
    } else {
      NA_real_
    },
    .groups = "drop"
  ) %>%
  filter(!is.na(P_Value), P_Value < 0.05)

tm_shape(p_value_stats, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 50), crs = 4326)) +
  tm_polygons(
    fill = "P_Value", 
    fill.scale = tm_scale(
      values = c("#00441b", "#006d2c", "#238b45", "#41ab5d", "#74c476"),
      limits = c(0, 0.05)
    ),
    fill.legend = tm_legend(title = "P-Value (< 0.05)", orientation = "portrait", position = tm_pos_in("right", "center")), 
    col = "white", 
    lwd = 0.05
  ) +
  tm_shape(world) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_compass(type = "arrow", size = 3, fontsize = 1, color.dark = "black", position = c("right", "top")) +
  tm_scalebar(text.size = 0.8, lwd = 3, color.dark = "black", position = c("left", "bottom")) +
  tm_title(text = "Statistically Significant Spatial Correlation", size = 1.5, fontface = "bold", position = c("center", "top")) +
  tm_layout(
    bg.color = "#e0f3f8", 
    legend.bg.color = "white",
    inner.margins = c(0.1, 0.05, 0.05, 0.05)
  )

## Scale bar set for latitude km and will be different at the top and bottom of the map.

grid_stats_temp <- joined_data %>%
  filter(!is.na(SST_Celsius), !is.na(Shear_kts)) %>%
  group_by(grid_id) %>%
  summarise(
    Count = n(),
    Correlation = if(n() > 5 && sd(SST_Celsius, na.rm=TRUE) > 0 && sd(Shear_kts, na.rm=TRUE) > 0) {
      cor(SST_Celsius, Shear_kts, use="complete.obs")
    } else {
      NA_real_
    },
    P_Value = if(n() > 5 && sd(SST_Celsius, na.rm=TRUE) > 0 && sd(Shear_kts, na.rm=TRUE) > 0) {
      cor.test(SST_Celsius, Shear_kts, method = "pearson")$p.value
    } else {
      NA_real_
    },
    .groups = "drop"
  ) %>%
  filter(!is.na(P_Value), P_Value < 0.05)

grid_stats <- grid_stats_temp %>%
  left_join(basin_lookup, by = "grid_id")

basin_stats <- grid_stats %>%
  st_drop_geometry() %>%
  group_by(Basin_Region) %>%
  summarise(
    Avg_r = round(mean(Correlation, na.rm = TRUE), 3), 
    Cells = n()
  ) %>%
  rename(Basin = Basin_Region)

tm_shape(grid_stats, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 50), crs = 4326)) +
  tm_polygons(
    fill = "Correlation", 
    fill.scale = tm_scale(
      values = c("#2c7bb6", "white", "#d7191c"), 
      midpoint = 0
    ),
    fill.legend = tm_legend(title = "Correlation (r)", orientation = "portrait", position = tm_pos_in("right", "center")), 
    col = NA,
    lwd = 0,
    fill_alpha = 0.6
  ) +
  tm_shape(world) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_compass(type = "arrow", size = 3, fontsize = 1, color.dark = "black", position = c("right", "top")) +
  tm_scalebar(text.size = 0.8, lwd = 3, color.dark = "black", position = c("left", "bottom")) +
  tm_title(text = "Correlation (r) & Significance (P < 0.05)", size = 1.5, fontface = "bold", position = c("center", "top")) +
  tm_layout(
    bg.color = "#e0f3f8", 
    legend.bg.color = "white",
    inner.margins = c(0.1, 0.05, 0.05, 0.05)
  )

## Scale bar set for latitude km and will be different at the top and bottom of the map.

ri_correlation_combined <- grid_stats_temp %>%
  left_join(
    ri_probability %>% st_drop_geometry() %>% select(grid_id, Probability), 
    by = "grid_id"
  )

high_ri_zones <- ri_correlation_combined %>% 
  filter(!is.na(Probability), Probability > 20) %>%
  mutate(RI_Zone = "RI Probability > 20%")

tm_shape(grid_stats, bbox = st_bbox(c(xmin = -100, xmax = -10, ymin = 5, ymax = 50), crs = 4326)) +
  tm_polygons(
    fill = "Correlation", 
    fill.scale = tm_scale(
      values = c("#2c7bb6", "white", "#d7191c"), 
      midpoint = 0
    ),
    fill.legend = tm_legend(title = "Correlation (r)", orientation = "portrait", position = tm_pos_in("right", "center")), 
    col = NA,
    lwd = 0,
    fill_alpha = 0.7
  ) +
  tm_shape(world) +
  tm_polygons(fill = "gray90", col = "gray60", lwd = 0.5) +
  tm_shape(high_ri_zones) +
  tm_borders(
    col = "RI_Zone",
    col.scale = tm_scale(values = "black"),
    col.legend = tm_legend(title = "", labels = "RI Probability > 20%", position = tm_pos_in("right", "center")),
    lwd = 1.2, 
    lty = "dashed"
  ) +
  tm_compass(type = "arrow", size = 3, fontsize = 1, color.dark = "black", position = c("right", "top")) +
  tm_scalebar(text.size = 0.8, lwd = 3, color.dark = "black", position = c("left", "bottom")) +
  tm_title(text = "Significant SST-Shear in High RI Regions", size = 1.5, fontface = "bold", position = c("center", "top")) +
  tm_layout(
    bg.color = "#e0f3f8", 
    legend.bg.color = "white",
    legend.frame = TRUE,
    inner.margins = c(0.1, 0.05, 0.05, 0.05)
  )

## The visual variable `col` of the layer "borders" and "polygons" contains a
## unique value. Therefore a categorical scale is applied (tm_scale_categorical).
## Scale bar set for latitude km and will be different at the top and bottom of the map.

ri_corr_summary <- ri_correlation_combined %>%
  st_drop_geometry() %>%
  filter(!is.na(Probability)) %>%
  mutate(RI_Category = case_when(
    Probability >= 30 ~ "Very High RI (≥30%)",
    Probability >= 20 ~ "High RI (20-30%)",
    Probability >= 10 ~ "Moderate RI (10-20%)",
    TRUE ~ "Low RI (<10%)"
  )) %>%
  group_by(RI_Category) %>%
  summarise(
    Avg_Correlation = round(mean(abs(Correlation), na.rm = TRUE), 3),
    Num_Cells = n(),
    Avg_Probability = round(mean(Probability, na.rm = TRUE), 1)
  ) %>%
  arrange(desc(Avg_Probability))

Question 6:

Congratulations! You just made spatial maps for your final project!!! Now that you have these maps, write 1-2 paragraphs about any trends or features that stand out to you on your maps from Question 4 and Question 5. Be detailed in your explanations. (15)

thx ;’’’’’’)

I think there are some really interesting, but not surprising, correlative relationships that are being displayed in the final geo-computed plots. To start, the incredibly low P-values displayed across the spatial domain confirms Kaplan and DeMaria (2003)’s analysis and thermodynamic TC theory that shear and SST are significantly correlated. As far as how far the correlative relationships go, the analysis in question 5 shows an interesting pattern of both positively correlated AND negatively correlated grids that are still statistically significant. This is an interesting trend that may be related to the graph of SST and Shear correlations found in our previous EDA assignment and proposal document that showed an awkward trend for mediocre shear and SST values to be positively correlated but strong or weak shear and SST values to be negatively correlated. As far as a deeper investigation into the cases of rapid intensification, thermodynamic TC theory is confirmed by the highest probabilities (RI%>20) occur on average in statistically significant grid spaces with with negative (mostly in the range from -1.0-0.6 r values, indicating extremely negative correlative relationships). This certainly makes sense intuitively; Rapid Intensification only occurs in environments where TCs have extremely favorable conditions (warm SSTs and low shear), so to see this ploted as result is not shocking.

Question 7:

Export your maps and Knit your HTML File. (5)

AI used (error in kniting):

“I tried knitting this file to an html but this was the render error I got.

processing file: Lab07_SpatialMapsPII.Rmd |………………………………………. | 89% [unnamed-chunk-4]

Error: ! object ‘grid_centroids’ not found Backtrace: ▆ 1. └─sf::st_coordinates(grid_centroids)

Quitting from Lab07_SpatialMapsPII.Rmd:421-627 [unnamed-chunk-4] Execution halted

What do I need to fix?”

Output:

“Looking at your code in Question 5, the error occurs because you’re trying to use grid_centroids on this line:

grid_coords <- st_coordinates(grid_centroids) But you never created grid_centroids! You need to add this line BEFORE you try to get the coordinates”