Lab5_CarissaFeliciano.qmd

---
title: "Lab 5"
author: "Carissa Feliciano"
format: html
embed-resources: true
---

# Prepare the data 
```{r}
library(data.table)
library(dtplyr)
library(dplyr)
library(DT)

#Load the data
if (!file.exists(file.path("~", "Downloads", "met_all.gz")))
  download.file(
    url = "https://raw.githubusercontent.com/USCbiostats/data-science-data/master/02_met/met_all.gz",
    destfile = file.path("~", "Downloads", "met_all.gz"),
    method   = "libcurl",
    timeout  = 60
    )
met <- fread(file.path("~", "Downloads", "met_all.gz"))
met <- as.data.frame(met)

if (!file.exists(file.path("~", "Downloads", "stations.gz")))
  download.file(
    url = "https://noaa-isd-pds.s3.amazonaws.com/isd-history.csv",
    destfile = file.path("~", "Downloads", "stations.gz"),
    method   = "libcurl",
    timeout  = 60
    )
stations <- fread(file.path("~", "Downloads", "stations.gz"))
stations <- as.data.frame(stations)
stations$USAF <- as.integer(stations$USAF)
```

```{r}
# Dealing with NAs and 999999s
stations$USAF[stations$USAF == 999999] <- NA
stations$CTRY[stations$CTRY == ""] <- NA
stations$STATE[stations$STATE == ""] <- NA
```

```{r}
# Selecting the three relevant columns and keeping unique records
stations <- unique(stations[, c('USAF', 'CTRY', 'STATE')])
```

```{r}
# Dropping the NAs
stations <- stations[!is.na('USAF'), ]
```

```{r}
# Removing duplicates
stations <- stations[!duplicated(stations$USAF), ]
```

```{r}
# Merge the data 
met <- merge(
  x     = met,      
  y     = stations, 
  by.x  = "USAFID",
  by.y  = "USAF", 
  all.x = TRUE,      
  all.y = FALSE
  )
```

# 1. Representative station for the US
```{r}
# Calculate medians 
median_temp <- quantile(met$temp, probs = 0.5, na.rm = TRUE)
median_wind <- quantile(met$wind.sp, probs = 0.5, na.rm = TRUE)
median_atm <- quantile(met$atm.press, probs = 0.5, na.rm = TRUE)

median_temp
median_wind
median_atm
```

```{r}
# Determine the three weather stations that best represent the continental US
met$dist_temp <- abs(met$temp - median_temp)
met$dist_wind <- abs(met$wind.sp - median_wind)
met$dist_atm <- abs(met$atm.press - median_atm)

met$dist_all <- met$dist_temp + met$dist_wind + met$dist_atm 

met[order(met$dist_all), ][1:3, ]
```

The three weather stations that best represent the continental US are stations 722246 (FL), 722390 (LA), and 723066 (NC). These three stations are all located in different states. 

```{r}
met <- met |> select (-dist_temp, -dist_wind, -dist_atm, -dist_all)
```

#2. Identify the most representative station per state.
```{r}
# Calculate the median for each variable by state
medians <- met |>
  group_by(STATE) |>
  summarize(
    med_temp = median(met$temp, na.rm = TRUE),
    med_wind = median(met$wind.sp, na.rm = TRUE),
    med_atm = median(met$atm.press, na.rm = TRUE)
  )

#Add the medians to the dataset 
met_medians <- met |>
  left_join(medians, by = "STATE")

# Calculate the euclidean distance from the median values
met_medians <- met_medians |>
  mutate(
    euc = sqrt(
      (temp - med_temp)^2 + 
      (wind.sp - med_wind)^2 +
      (atm.press - med_atm)^2
    )
  )
```

```{r}
# Find the stations closest to the median per state
median_stations <- met_medians |>
  group_by(STATE) |>
  filter(euc == min(euc, na.rm = TRUE)) |>
  slice(which.min(lat)) |>
  ungroup() |>
  select(USAFID, STATE, temp, wind.sp, atm.press, med_temp, med_wind, med_atm, lat, lon)

nrow(median_stations)

datatable(median_stations, options = list(scrollX = TRUE, pageLength = 10))
```
The representative stations per state are displayed in the table above. Representative stations are missing for ND and WA because there are no stations in these states that have values recorded for temperature, wind speed, and atmospheric pressure.

# 3. Identify the station closest to the mid-point of the state. 
Determine the mid-point for each state.
```{r}
library(dplyr)
library(ggplot2)

# Get the map data for US states
states <- map_data("state")

# Calculate the mid-point for each state
midpoints <- states |>
  group_by(region) |>
  summarize(
    mid_lon = mean(long, na.rm = TRUE),
    mid_lat = mean(lat, na.rm = TRUE)
  )

# Replace region column with lowercase state names with abbreviations
state_midpoints <- midpoints |> 
  mutate(STATE = state.abb[match(region, tolower(state.name))])

state_midpoints <- state_midpoints |>
  select(-region)
  
# Merge midpoint data with met dataset 
met_midpoints <- merge(met, state_midpoints, by = "STATE")
```

Identify the station closest to the mid-point for each state. 
```{r}
# Calculate the distance from the stations to the mid-point
library(geosphere)

met_midpoints <- met_midpoints |>
  mutate(dist = distHaversine(cbind(lon,lat), cbind(mid_lon, mid_lat)))
```

```{r}
midpoint_stations <- met_midpoints |>
  group_by(STATE) |>
  filter(dist == min(dist, na.rm = TRUE)) |>
  ungroup() |>
  distinct(USAFID, .keep_all = TRUE) |>
  select(STATE, USAFID, lat, lon, mid_lon, mid_lat, dist)

nrow(midpoint_stations)

datatable(midpoint_stations, options = list(scrollX = TRUE, pageLength = 10))
```

Create a leaflet map of all points
```{r}
# Add label column to median stations dataset
median_stations$category <- "Median"

# Add label column to midpoint stations dataset 
midpoint_stations$category <- "Mid-point"

# Merge columns 
all_stations <- bind_rows(
  median_stations |> select(STATE, USAFID, lat, lon, category),
  midpoint_stations |> select(STATE, USAFID, lat, lon, category)
)

nrow(all_stations)
```

```{r}
#Plot leaflet map
library(leaflet)
pal <- colorFactor(c("blue", "orange"), domain = c("Median", "Mid-point"))

leaflet(all_stations) |>
  addTiles() |>
  addCircleMarkers(
    lng = ~lon, 
    lat = ~lat, 
    color = ~pal(category),
    stroke = FALSE,
    opacity = 0.7, fillOpacity = 0.7, radius = 5,
    popup = ~paste("Station ID:", USAFID, "<br> Category:", category)
  ) |>
  addLegend("bottomleft",
            pal = pal,
            values = ~category,
            title = "Category",
            opacity = 1)
  
```

# 4. Means of means
```{r}
# Compute the states' average temperature 
met <- met |>
  group_by(STATE) |>
  mutate(
    mean_temp = mean(temp, na.rm = TRUE)) |>
  ungroup()
  
# Use the state avg temp to classify them 
met$temp_cat <- rep(NA, nrow(met))
met$temp_cat[met$mean_temp <20] <- 'low'
met$temp_cat[met$mean_temp >= 20 & met$mean_temp <25] <- 'mid'
met$temp_cat[met$mean_temp >= 25] <- 'high'
```
```{r}
summary_table <- met |>
  group_by(temp_cat) |>
  summarize(
    `Number of entries` = n(),
    `Number of NA entries` = sum(is.na(temp) | is.na(wind.sp) | is.na(atm.press)),
    `Number of stations` = n_distinct(USAFID),
    `Number of states` = n_distinct(STATE),
    `Mean temp` = mean(temp, na.rm = TRUE),
    `Mean wind speed` = mean(wind.sp, na.rm = TRUE),
    `Mean atm pressure` = mean(atm.press, na.rm = TRUE)
  ) |>
  mutate(temp_cat = factor(temp_cat, levels = c("low", "mid", "high"))) |>
  arrange(temp_cat)
    
datatable(summary_table, options = list(scrollX = TRUE, pageLength = 3))
```