Skip to content

Commit

Permalink
Work on calc_outbreak_history. Remove circle geometry and define max …
Browse files Browse the repository at this point in the history
…distance allowed in distance matrix.
  • Loading branch information
n8layman committed Jun 27, 2024
1 parent f14b830 commit 76b2670
Showing 1 changed file with 31 additions and 23 deletions.
54 changes: 31 additions & 23 deletions R/calc_outbreak_history.R
Original file line number Diff line number Diff line change
Expand Up @@ -44,21 +44,21 @@ calc_outbreak_history <- function(wahis_rvf_outbreaks_preprocessed,
recent = 1/6) {

recent_outbreaks <- wahis_rvf_outbreaks_preprocessed |>
mutate(end_date = coalesce(min(current_date, outbreak_end_date, na.rm = T), outbreak_start_date),
mutate(end_date = pmin(current_date, coalesce(outbreak_end_date, outbreak_start_date), na.rm = T),
years_since = as.numeric(as.duration(current_date - end_date), "years"),
weight = ifelse(is.na(cases) | cases == 1, 1,log10(cases + 1))*exp(-beta_time*years_since)) |>
filter(years_since < recent & years_since > 0) |>
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

old_outbreaks <- wahis_rvf_outbreaks_preprocessed |>
mutate(end_date = coalesce(min(current_date, outbreak_end_date, na.rm = T), outbreak_start_date),
mutate(end_date = pmin(current_date, coalesce(outbreak_end_date, outbreak_start_date), na.rm = T),
years_since = as.numeric(as.duration(current_date - end_date), "years"),
weight = ifelse(is.na(cases) | cases == 1, 1,log10(cases + 1))*exp(-beta_time*years_since)) |>
filter(years_since >= recent & years_since < max_years) |>
sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)

raster_template <- terra::rast(continent_raster_template)
raster_template[] <- 0
raster_template[] <- 1
raster_template <- crop(raster_template, vect(continent_polygon$geometry), mask = TRUE)

recent_outbreak_weights <- get_outbreak_weights(recent_outbreaks, raster_template, )
Expand All @@ -79,27 +79,23 @@ calc_outbreak_history <- function(wahis_rvf_outbreaks_preprocessed,
# of tibble with lat, long, recent_outbreak_weight, old_outbreak_weight
}


get_outbreak_weights <- function(outbreaks, raster_template, ) {

if(nrow(outbreaks) == 0) {
raster[] <- 0
return(raster)
}
# This needs to happen for each outbreak separately otherwise when we calculate weights
# They can blend together. What I mean is that the distance matrix contains distances between
# cells within each circle and every origin not just the origin of that circle.
get_outbreak_weights <- function(outbreaks, points, origin, ) {

outbreak_circles <- vect(outbreaks$geometry) |>
terra::buffer(within_km * 1000)
# outbreak_circles <- vect(outbreaks$geometry) |>
# terra::buffer(within_km * 1000)
# Unpack the raster that we want to calculate

# raster <- crop(raster_template, outbreak_circles, mask = TRUE)
raster <- terra::mask(raster_template, outbreak_circles)
# raster <- terra::mask(raster_template, outbreak_circles, updatevalue = 0)
# For each pixel in the raster, calculate the most recent outbreak within `within_km` km
xy <- s2::as_s2_lnglat(terra::crds(raster))
idx <- which(!is.nan(raster[]))
xy <- s2::as_s2_lnglat(terra::crds(raster_template))
idx <- which(!is.nan(raster_template[]))

#Get the unique outbreaks$geometry values while maintaining the sf object characteristics


# Create an index for the unique values of outbreaks$geometry and a column of that index for each outbreak
locations <- unique(outbreaks$geometry)
attributes(locations) <- attributes(outbreaks$geometry)
Expand All @@ -108,14 +104,26 @@ get_outbreak_weights <- function(outbreaks, raster_template, ) {
xy_o <- s2::as_s2_lnglat(locations)

# For each pixel identify the outbreaks within `within_km` km
#matches <- s2::s2_dwithin_matrix(xy, xy_o, within_km * 1000)
# matches <- s2::s2_dwithin_matrix(xy, xy_o, within_km * 1000)
# 883 origins (xy_o) x 133491 cells (xy) = 117872553 distances.
# but some of those distances cross the buffer. It's every cell to every origin.
# We don't need to calculate that just within the circles around each origin
# Also this doesn't connect to log cases at all. It's just a strict spatial
# distance weight at this point.
matches_dist <- s2::s2_distance_matrix(xy, xy_o)
weights <- rowSums(exp(-beta_dist*matches_dist/1000))
raster[idx] <- weights
test <- terra::cover(raster, raster_template)

as.data.frame(test, xy = T)

# Convert raster into lat-long tibble. Needs all pixels
# Drop all distances greater than within_km
matches_dist[matches_dist > within_km] <- NA

# Calculate a weighting factor based on distance. Note we haven't included log10 cases yet.
# This is negative exponential decay - points closer to the origin will be 1 and those farther
# away will be closer to zero mediated by beta_dist.
# Not sure about this rowSums bit yet.
weights <- rowSums(exp(-beta_dist*matches_dist/1000)) # /1000 to get to meters

raster_template[idx] <- weights

# Convert xy into tibble and add weights column
# as_tibble(xy) |> rename(long = x, lat = y) |> mutate(weight = weights)

}

0 comments on commit 76b2670

Please sign in to comment.