R/epitool_spatial.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# =============================================================================
# R4VN EpiTool Studio - Spatial epidemiology engines
# File: R/epitool_spatial.R
# =============================================================================

.r4vn_epi_require_sf <- function() {
  if(!requireNamespace("sf",quietly=TRUE)) {
    stop("Spatial epidemiology requires the optional package `sf`.",call.=FALSE)
  }
  invisible(TRUE)
}

.r4vn_epi_utm_epsg <- function(lon,lat) {
  lon <- mean(lon[is.finite(lon)],na.rm=TRUE); lat <- mean(lat[is.finite(lat)],na.rm=TRUE)
  if(!is.finite(lon)||!is.finite(lat)) return(3857L)
  zone <- floor((lon+180)/6)+1
  zone <- max(1,min(60,zone))
  if(lat>=0) as.integer(32600+zone) else as.integer(32700+zone)
}

.r4vn_epi_spatial_detect <- function(data) {
  p <- .r4vn_epi_profile(data)
  list(
    latitude=if(length(p$lat_candidates))p$lat_candidates[1] else NA_character_,
    longitude=if(length(p$lon_candidates))p$lon_candidates[1] else NA_character_,
    id=if(length(p$id_candidates))p$id_candidates[1] else NA_character_,
    onset=if(length(p$dates))p$dates[1] else NA_character_,
    status=if(length(p$case_candidates))p$case_candidates[1] else NA_character_
  )
}

.r4vn_epi_spatial_validate <- function(data,lat,lon) {
  if(!all(c(lat,lon)%in%names(data))) stop("Latitude or longitude variable was not found.",call.=FALSE)
  la <- suppressWarnings(as.numeric(data[[lat]]))
  lo <- suppressWarnings(as.numeric(data[[lon]]))
  valid <- is.finite(la)&is.finite(lo)&la>=-90&la<=90&lo>=-180&lo<=180
  swapped <- is.finite(la)&is.finite(lo)&abs(la)>90&abs(lo)<=90
  duplicate_xy <- duplicated(data.frame(la,lo)) & valid
  list(
    latitude=la,longitude=lo,valid=valid,
    n=nrow(data),valid_n=sum(valid),missing_n=sum(!is.finite(la)|!is.finite(lo)),
    invalid_range_n=sum(is.finite(la)&is.finite(lo)&!valid),
    possible_swapped_n=sum(swapped),
    duplicate_location_n=sum(duplicate_xy),
    completeness=if(nrow(data)>0)sum(valid)/nrow(data) else NA_real_
  )
}

.r4vn_epi_as_sf <- function(data,lat,lon,remove=FALSE,crs=4326) {
  .r4vn_epi_require_sf()
  v <- .r4vn_epi_spatial_validate(data,lat,lon)
  d <- data[v$valid,,drop=FALSE]
  sf::st_as_sf(d,coords=c(lon,lat),crs=crs,remove=remove)
}

.r4vn_epi_spatial_distance <- function(cases,sources,
                                       case_lat="latitude",case_lon="longitude",
                                       source_lat="latitude",source_lon="longitude",
                                       source_id="source_id") {
  .r4vn_epi_require_sf()
  c_sf <- .r4vn_epi_as_sf(cases,case_lat,case_lon,remove=FALSE)
  s_sf <- .r4vn_epi_as_sf(sources,source_lat,source_lon,remove=FALSE)
  if(!nrow(c_sf)||!nrow(s_sf)) stop("Valid case and source coordinates are required.",call.=FALSE)
  ctr <- sf::st_coordinates(sf::st_centroid(sf::st_union(sf::st_geometry(c_sf))))
  epsg <- .r4vn_epi_utm_epsg(ctr[1],ctr[2])
  cp <- sf::st_transform(c_sf,epsg); sp <- sf::st_transform(s_sf,epsg)
  idx <- sf::st_nearest_feature(cp,sp)
  d <- as.numeric(sf::st_distance(cp,sp[idx,],by_element=TRUE))
  id <- if(source_id%in%names(sources)) as.character(sources[[source_id]]) else as.character(seq_len(nrow(sources)))
  out <- sf::st_drop_geometry(c_sf)
  out$nearest_source <- id[idx]
  out$distance_source_m <- d
  out
}

.r4vn_epi_spatial_buffer <- function(cases,sources,radius_m=500,
                                     case_lat="latitude",case_lon="longitude",
                                     source_lat="latitude",source_lon="longitude",
                                     source_id="source_id") {
  .r4vn_epi_require_sf()
  if(!is.finite(radius_m)||radius_m<=0) stop("`radius_m` must be positive.",call.=FALSE)
  c_sf <- .r4vn_epi_as_sf(cases,case_lat,case_lon,remove=FALSE)
  s_sf <- .r4vn_epi_as_sf(sources,source_lat,source_lon,remove=FALSE)
  allpts <- c(sf::st_coordinates(c_sf)[,1],sf::st_coordinates(s_sf)[,1])
  allats <- c(sf::st_coordinates(c_sf)[,2],sf::st_coordinates(s_sf)[,2])
  epsg <- .r4vn_epi_utm_epsg(allpts,allats)
  cp <- sf::st_transform(c_sf,epsg); sp <- sf::st_transform(s_sf,epsg)
  buffers <- sf::st_buffer(sp,dist=radius_m)
  inside_list <- sf::st_intersects(cp,buffers)
  nearest <- sf::st_nearest_feature(cp,sp)
  dist_nearest <- as.numeric(sf::st_distance(cp,sp[nearest,],by_element=TRUE))
  source_labels <- if(source_id%in%names(sources)) as.character(sources[[source_id]]) else as.character(seq_len(nrow(sources)))
  out <- sf::st_drop_geometry(c_sf)
  out$within_buffer <- lengths(inside_list)>0
  out$n_sources_within_buffer <- lengths(inside_list)
  out$nearest_source <- source_labels[nearest]
  out$distance_source_m <- dist_nearest
  out$buffer_radius_m <- radius_m
  list(
    data=out,
    cases_sf=c_sf,
    sources_sf=s_sf,
    buffers_sf=sf::st_transform(buffers,4326),
    radius_m=radius_m,
    projected_crs=epsg
  )
}

.r4vn_epi_spatial_multi_buffer <- function(cases,sources,radii_m=c(250,500,1000,2000),
                                           case_lat="latitude",case_lon="longitude",
                                           source_lat="latitude",source_lon="longitude",
                                           source_id="source_id") {
  radii_m <- sort(unique(as.numeric(radii_m[is.finite(radii_m)&radii_m>0])))
  if(!length(radii_m)) stop("Provide at least one positive radius.",call.=FALSE)
  z <- .r4vn_epi_spatial_distance(cases,sources,case_lat,case_lon,source_lat,source_lon,source_id)
  br <- c(-Inf,radii_m,Inf)
  labels <- c(paste0("\u2264",radii_m[1]," m"),
              if(length(radii_m)>1) paste0(radii_m[-length(radii_m)]+1,"\u2013",radii_m[-1]," m") else character(),
              paste0(">",tail(radii_m,1)," m"))
  z$distance_band <- cut(z$distance_source_m,breaks=br,labels=labels,include.lowest=TRUE,right=TRUE)
  z
}

.r4vn_epi_spatial_grid <- function(data,lat,lon,cell_m=500,shape=c("hexagon","square")) {
  .r4vn_epi_require_sf()
  shape <- match.arg(shape)
  pts <- .r4vn_epi_as_sf(data,lat,lon,remove=FALSE)
  if(!nrow(pts)) stop("No valid coordinates.",call.=FALSE)
  xy <- sf::st_coordinates(pts)
  epsg <- .r4vn_epi_utm_epsg(xy[,1],xy[,2])
  pp <- sf::st_transform(pts,epsg)
  grid <- sf::st_make_grid(pp,cellsize=cell_m,square=(shape=="square"))
  grid <- sf::st_sf(grid_id=seq_along(grid),geometry=grid)
  rel <- sf::st_intersects(grid,pp)
  grid$cases <- lengths(rel)
  grid <- grid[grid$cases>0,,drop=FALSE]
  sf::st_transform(grid,4326)
}

.r4vn_epi_spatial_dbscan <- function(data,lat,lon,eps_m=500,minPts=5L) {
  .r4vn_epi_require_sf()
  if(!requireNamespace("dbscan",quietly=TRUE)) stop("Point clustering requires the optional package `dbscan`.",call.=FALSE)
  pts <- .r4vn_epi_as_sf(data,lat,lon,remove=FALSE)
  xy0 <- sf::st_coordinates(pts)
  epsg <- .r4vn_epi_utm_epsg(xy0[,1],xy0[,2])
  pp <- sf::st_transform(pts,epsg)
  xy <- sf::st_coordinates(pp)
  fit <- dbscan::dbscan(xy,eps=eps_m,minPts=as.integer(minPts))
  out <- sf::st_drop_geometry(pts)
  out$cluster_id <- fit$cluster
  out$is_clustered <- fit$cluster>0
  attr(out,"dbscan") <- fit
  out
}

.r4vn_epi_spatial_join_area <- function(data,boundary,lat,lon,fields=NULL) {
  .r4vn_epi_require_sf()
  if(!inherits(boundary,"sf")) stop("`boundary` must be an sf object.",call.=FALSE)
  pts <- .r4vn_epi_as_sf(data,lat,lon,remove=FALSE)
  if(is.na(sf::st_crs(boundary))) stop("Boundary file has no coordinate reference system.",call.=FALSE)
  pts <- sf::st_transform(pts,sf::st_crs(boundary))
  if(is.null(fields)) fields <- setdiff(names(boundary),attr(boundary,"sf_column") %||% "geometry")
  fields <- intersect(fields,names(boundary))
  joined <- sf::st_join(pts,boundary[,fields,drop=FALSE],left=TRUE,join=sf::st_within)
  sf::st_drop_geometry(joined)
}

`%||%` <- function(x,y) if(is.null(x)||!length(x)) y else x

.r4vn_epi_spatial_area_rates <- function(cases,boundary,area_boundary,
                                         population,area_population,
                                         multiplier=100000,
                                         case_lat=NULL,case_lon=NULL) {
  .r4vn_epi_require_sf()
  if(!inherits(boundary,"sf")) stop("Boundary must be an sf object.",call.=FALSE)
  if(!area_boundary%in%names(boundary)) stop("Boundary area ID not found.",call.=FALSE)
  if(!area_population%in%names(population)) stop("Population area ID not found.",call.=FALSE)
  if(!"population"%in%names(population)) stop("Population data must contain a `population` field after mapping.",call.=FALSE)

  if(inherits(cases,"sf")) {
    pts <- sf::st_transform(cases,sf::st_crs(boundary))
    joined <- sf::st_join(pts,boundary[,area_boundary,drop=FALSE],left=TRUE,join=sf::st_within)
    cc <- as.data.frame(table(sf::st_drop_geometry(joined)[[area_boundary]]),stringsAsFactors=FALSE)
    names(cc) <- c("area_id","cases")
  } else if(!is.null(case_lat)&&!is.null(case_lon)) {
    pts <- .r4vn_epi_as_sf(cases,case_lat,case_lon,remove=FALSE)
    pts <- sf::st_transform(pts,sf::st_crs(boundary))
    joined <- sf::st_join(pts,boundary[,area_boundary,drop=FALSE],left=TRUE,join=sf::st_within)
    cc <- as.data.frame(table(sf::st_drop_geometry(joined)[[area_boundary]]),stringsAsFactors=FALSE)
    names(cc) <- c("area_id","cases")
  } else if("cases"%in%names(cases)) {
    cc <- data.frame(area_id=as.character(cases[[area_boundary]]),cases=as.numeric(cases$cases))
  } else stop("Cases must be points with coordinates or area-level counts.",call.=FALSE)

  pop <- data.frame(area_id=as.character(population[[area_population]]),
                    population=suppressWarnings(as.numeric(population$population)),
                    stringsAsFactors=FALSE)
  tab <- merge(pop,cc,by="area_id",all.x=TRUE); tab$cases[is.na(tab$cases)] <- 0
  tab$rate <- ifelse(tab$population>0,tab$cases/tab$population*multiplier,NA_real_)
  b <- boundary
  b$.area_join <- as.character(b[[area_boundary]])
  out <- merge(b,tab,by.x=".area_join",by.y="area_id",all.x=TRUE,sort=FALSE)
  out$multiplier <- multiplier
  out
}

.r4vn_epi_spatial_weights <- function(polygons,style="W",queen=TRUE) {
  if(!requireNamespace("spdep",quietly=TRUE)) stop("Hotspot analysis requires the optional package `spdep`.",call.=FALSE)
  nb <- spdep::poly2nb(polygons,queen=queen)
  spdep::nb2listw(nb,style=style,zero.policy=TRUE)
}

.r4vn_epi_spatial_local_moran <- function(polygons,value,alpha=.05,queen=TRUE) {
  .r4vn_epi_require_sf()
  if(!value%in%names(polygons)) stop("Value variable not found.",call.=FALSE)
  x <- suppressWarnings(as.numeric(polygons[[value]]))
  if(any(!is.finite(x))) stop("Local Moran analysis requires finite numeric values for every area.",call.=FALSE)
  lw <- .r4vn_epi_spatial_weights(polygons,queen=queen)
  lm <- spdep::localmoran(x,lw,zero.policy=TRUE,na.action=na.exclude)
  lagx <- spdep::lag.listw(lw,x,zero.policy=TRUE)
  z <- as.numeric(scale(x)); lz <- as.numeric(scale(lagx))
  pcol <- grep("^Pr",colnames(lm),value=TRUE)[1]
  p <- if(length(pcol)) lm[,pcol] else rep(NA_real_,length(x))
  cls <- ifelse(p>alpha|!is.finite(p),"Not significant",
         ifelse(z>=0&lz>=0,"High-High",
         ifelse(z<0&lz<0,"Low-Low",
         ifelse(z>=0&lz<0,"High-Low","Low-High"))))
  out <- polygons
  out$local_moran_I <- lm[,1]
  out$local_moran_p <- p
  out$local_moran_cluster <- factor(cls,levels=c("High-High","Low-Low","High-Low","Low-High","Not significant"))
  out
}

.r4vn_epi_spatial_getis <- function(polygons,value,alpha=.05,queen=TRUE) {
  .r4vn_epi_require_sf()
  if(!value%in%names(polygons)) stop("Value variable not found.",call.=FALSE)
  x <- suppressWarnings(as.numeric(polygons[[value]]))
  if(any(!is.finite(x))) stop("Getis-Ord analysis requires finite numeric values for every area.",call.=FALSE)
  if(!requireNamespace("spdep",quietly=TRUE)) stop("Getis-Ord analysis requires the optional package `spdep`.",call.=FALSE)
  nb <- spdep::poly2nb(polygons,queen=queen)
  nb_star <- spdep::include.self(nb)
  lw <- spdep::nb2listw(nb_star,style="B",zero.policy=TRUE)
  g <- spdep::localG(x,lw,zero.policy=TRUE)
  z <- as.numeric(g)
  p <- 2*stats::pnorm(-abs(z))
  cls <- ifelse(p>alpha,"Not significant",ifelse(z>0,"Hotspot","Coldspot"))
  out <- polygons
  out$getis_g_z <- z
  out$getis_p <- p
  out$getis_class <- factor(cls,levels=c("Hotspot","Coldspot","Not significant"))
  out
}

.r4vn_epi_spatial_privacy <- function(data,lat,lon,mode=c("analysis","presentation"),
                                      jitter_m=150,round_digits=NULL,seed=NULL) {
  mode <- match.arg(mode)
  out <- data
  if(mode=="analysis") return(out)
  if(!is.null(round_digits)) {
    out[[lat]] <- round(as.numeric(out[[lat]]),round_digits)
    out[[lon]] <- round(as.numeric(out[[lon]]),round_digits)
  }
  if(is.finite(jitter_m)&&jitter_m>0) {
    .r4vn_set_seed_if(seed)
    la <- as.numeric(out[[lat]]); lo <- as.numeric(out[[lon]])
    # Approximate displacement for display only; engines must use original data.
    angle <- runif(length(la),0,2*pi); r <- sqrt(runif(length(la),0,1))*jitter_m
    dlat <- (r*cos(angle))/111320
    dlon <- (r*sin(angle))/(111320*pmax(.1,cos(la*pi/180)))
    out[[lat]] <- la+dlat; out[[lon]] <- lo+dlon
  }
  out
}

.r4vn_epi_spatial_summary <- function(data,lat,lon) {
  v <- .r4vn_epi_spatial_validate(data,lat,lon)
  data.frame(
    Metric=c("Records","Valid coordinates","Spatial completeness","Missing coordinates",
             "Invalid coordinates","Possible reversed coordinates","Duplicate locations"),
    Value=c(v$n,v$valid_n,.r4vn_epi_pct(v$completeness),v$missing_n,
            v$invalid_range_n,v$possible_swapped_n,v$duplicate_location_n),
    stringsAsFactors=FALSE
  )
}