In UW-Upwelling-Project/imageryML: Functions for Machine Learning with Oceanographic Imagery

knitr::opts_chunk$set(out.width = "70%", fig.align = "center")

This example shows how work through the points along the 20km coastline and return the upwelling for each day in the year. See vignette("auto3") for how to do this for one year. This vignette is just the code.

Load the needed packages for plotting and data.

require(raster)
require(imageryML) # the data
require(magrittr) # for %>%

Download the data

Set region lat-lon of interest.

min_lat <- 40.625
max_lat <- 50.625
min_lon <- 229.875
max_lon <- 236.625

If the file does not exist, then download.

for (yr in 1982:2020) {
  start_date <- paste0(yr, "-01-01")
  stop_date <- paste0(yr, "-12-31")
  file_name <- file.path(here::here(), "inst", "extdata", paste0(
    "erddap_", start_date, "_", stop_date, "_", "lat_", min_lat, "_", max_lat,
    "_lon_", min_lon, "_", max_lon, ".csv"
  ))
  if (file.exists(file_name)) next
  cat(yr, " ")
  import_erddap(start_date, stop_date,
    min_lat = min_lat, max_lat = max_lat,
    min_lon = min_lon, max_lon = max_lon,
    silent = TRUE
  )
}

Process raw data

Notice that I substract 360 off the longitudes. This is crucial. The standard projections use this location for the 0 longitude line.

processDF <- function(df_raw) {
  df_processed <- df_raw %>%
    # Get rid of miscellaneous zlev in first row
    dplyr::slice(-1) %>%
    # zlev is a column of zeroes, so get rid of that
    dplyr::select(-zlev) %>%
    # Convert into date
    dplyr::mutate(
      time = lubridate::ymd_hms(time),
      # Extract out day so I can just filter for first day in each month
      day = lubridate::day(time)
    ) %>%
    dplyr::select(-day) %>%
    # Set column names
    dplyr::rename(
      date = time,
      lat = latitude,
      lon = longitude
    ) %>%
    # Convert date column to Date type
    dplyr::mutate(
      date = as.Date(date),
      lat = as.numeric(lat),
      lon = as.numeric(lon) - 360,
      sst = as.numeric(sst)
    )
  return(df_processed)
}

Set up needed lines and raster

Need to create a sample raster on the Winkel Tripel projection (what the buffer 300km line is in) and then get its bounding box and crop our 300km line to that. vignette("data") describes the 300km, 20km and sample points data sets.

efil <- "erddap_2010-01-01_2010-12-31_lat_40.625_50.625_lon_229.875_236.625.csv"
filename = system.file("extdata", efil, package = "imageryML")
df_raw <- readr::read_csv(filename)
df_processed <- processDF(df_raw)
custom_date <- "2010-08-04"
x <- df_processed %>% dplyr::filter(date == custom_date)
r <- x %>% dplyr::select(lon, lat, sst) %>% raster::rasterFromXYZ(crs="+proj=longlat")

# Re-project the SST raster into Winkel Tripel projection (which is what 300km line is in)
wintri.crs <- raster::crs(sample_points[["wintri"]]$km100, asText = TRUE)
r.wintri <- raster::projectRaster(r, crs = wintri.crs, over = TRUE)

# Crop down the 300km line
r.wintri.bbox <- raster::bbox(r.wintri)
l.offshore <- smoothr::smooth(buffer300$wintri$line, method = "ksmooth", smoothness = 5)
l.offshore <- raster::crop(l.offshore, r.wintri.bbox)

Run through the years

final_df <- c()
for (yr in 1982:2020) {
  cat(yr, " ")
  start_date <- paste0(yr, "-01-01")
  stop_date <- paste0(yr, "-12-31")
  efil <- paste0("erddap_", start_date, "_", stop_date, "_", "lat_", min_lat, "_", max_lat, "_lon_", min_lon, "_", max_lon, ".csv")
  filename <- system.file("extdata", efil, package = "imageryML")
  df_raw <- readr::read_csv(filename, show_col_types = FALSE)
  df_processed <- processDF(df_raw)
  dates <- unique(df_processed$date)
  for (i in 1:length(dates)) {
    x <- df_processed %>% dplyr::filter(date == dates[i])
    out <- autoDetect3(x, l.offshore = l.offshore, threshold = 2, val = "sst", smooth.method = "none")$df
    out <- cbind(date = dates[i], out)
    final_df <- rbind(final_df, out)
  }
}
final_df <- as.data.frame(final_df)

Add on the lat-lon values

For making plots, it will be helpful to have the lat-lon values for the coast points.

# 1. make the x and y into spatial points. need to use the wintri projection
tmp <- sp::SpatialPoints(final_df %>% dplyr::select(x, y), proj4string = raster::crs(wintri.crs))
# 2. transform into the +proj=longlat; I am using crs(r) just to make sure I use the same crs as the original raster (which should be "+proj=longlat")
tmp <- sp::spTransform(tmp, raster::crs(r))
# 3. Add the lat-lon onto final_df
final_df_wlatlon <-cbind(final_df, lon=tmp@coords[,1], lat=tmp@coords[,2])
final_df_wlatlon <- final_df_wlatlon[c("date", "x", "y", "lon", "lat", 
                       "offshore.SST", "coast.SST", "upwelling")]

# Not for users. This is for adding the file to the extdata dir.
file_name <- file.path(here::here(), "inst", "extdata", "final_df_wlatlon.rda")
save(final_df_wlatlon, file=file_name)

Summary

final_df_wlatlon now has the upwelling logical value for each day from 1982 to 2020. This data frame is save in extdata as final_df_wlatlon.rda. This data.frame can be used for making heatmaps and contour maps of the upwelling timing.

Available at

filename = system.file("extdata", "final_df_wlatlon.rda", package = "imageryML")

UW-Upwelling-Project/imageryML documentation built on Dec. 18, 2021, 6:11 p.m.