R/simplify.R

In heike/ggmapr: Manipulation of Maps

#' Douglas Pecker algorithm for line thinning
#'
#' Implementation of the Douglas-Peucker algorithm for line thinning
#' http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
#' @param x vector of longitudes
#' @param y vector of latitudes
#' @param tol tolerance
simplify_rec <- function(x, y, tol = 0.01) {
  n <- length(x)
  if (n <= 2) return(unique(c(1, n)))
  dist <- point_line_dist(x, y, x[1], y[1], x[n], y[n])

  if (max(dist, na.rm = T) > tol) {
    furthest <- which.max(dist)
    c(
      simplify_rec(x[1:(furthest)], y[1:(furthest)], tol),
      # don't include furthest into list of essential points twice, but base the computation on it
      simplify_rec(x[(furthest):n], y[(furthest):n], tol)[-1] + furthest-1
    )
  } else {
    return(c(1,n))
  }
}

simplify_poly <- function(x, y, tol = 0.01) {
  n <- length(x)
  if (n <= 2) return(unique(c(1, n)))
  dist <- point_line_dist(x, y, x[1], y[1], x[n %/% 2], y[n %/% 2])
  if (max(dist, na.rm = T) > tol) {
    furthest <- which.max(dist)
    c(
      simplify_rec(x[1:(furthest)], y[1:(furthest)], tol),
      # don't include furthest into list of essential points twice, but base the computation on it
      simplify_rec(x[(furthest):n], y[(furthest):n] , tol)[-1] + furthest-1
    )
  } else {
    return(c(1,n))
  }
}

#' Distance between point and line
#'
#' Compute distance between point given as (px, py) and line spanned by points (lx1, ly1) and (lx2, ly2).
#' From http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
#' @param px x coordinate of  point outside the
#' @param py y coordinate of  point
#' @param lx_1, x  coordinate of 1st point spanning a line
#' @param ly_1, y  coordinate of 1st point spanning a line
#' @param lx_2, x  coordinate of 2nd point spanning a line
#' @param ly_2, y  coordinate of 2nd point spanning a line
point_line_dist <- function(px, py, lx_1, ly_1, lx_2, ly_2) {
  abs((lx_2 - lx_1) * (ly_1 - py) - (lx_1 - px) * (ly_2 - ly_1)) /
    sqrt((lx_2 - lx_1) ^ 2 + (ly_2 - ly_1) ^ 2)
}

# Precompute all tolerances so that we can post-select quickly
compute_tol <- function(points, offset = 0) {
  n <- nrow(points)
  if (n <= 2) {
    c()
  } else if (n == 3) {
    with(points,
         point_line_dist(long[2], lat[2], long[1], lat[1], long[3], lat[3]))
  } else {
    dist <- with(points,
                 point_line_dist(long[2:(n-1)], lat[2:(n-1)], long[1], lat[1], long[n], lat[n])
    )

    furthest <- which.max(dist)
    if (length(furthest) == 0) browser()
    c(
      compute_tol(points[1:(furthest + 1), ], offset),
      dist[furthest],
      compute_tol(points[(furthest + 1):n, ], furthest + offset)
    )
  }
}

#' Simplify a polygon shape by straightening out jittery lines
#'
#' @param map map object as generated by functions such as `process_shape`, `ggplot2::map_data`, or `ggplot2::fortify`
#' @param tol tolerance to use for simplifying the map. Larger tolerances will lead to smaller maps at lower resolutions.
#' @export
#' @examples
#' states05 <- simplify(states, 0.001)
#' states05 %>% ggplot(aes(x = long, y = lat)) + geom_path(aes(group = group))
#' states01 <- simplify(states, 0.01)
#' states01 %>% ggplot(aes(x = long, y = lat)) + geom_path(aes(group = group))
simplify <- function(map, tol = 0.01) {
  group <- data <- NULL

  mlist <- map %>% tidyr::nest(-group)
  mlist$data <- mlist$data %>% purrr::map(
    .f = function(d) {
      res <- simplify_poly(x = d$long, y = d$lat, tol = tol)
      d[res,]
    }
  )
  tidyr::unnest(mlist)
}