R/sf_ext.R

In faridcher/futils: Farid utility functions

#' @export
plot.sfc_LINESTRING <- function(
  x,
  y,
  ...,
  cex = 1,
  lty = 1,
  lwd = 1,
  col = 1,
  pch = 1,
  type = 'l',
  add = FALSE) {
  # I added cex param
  # FIXME: take care of lend, ljoin, and lmitre
  stopifnot(missing(y))
  if (! add)
    plot_sf(x, ...)
  lty = rep(lty, length.out = length(x))
  lwd = rep(lwd, length.out = length(x))
  col = rep(col, length.out = length(x))
  pch  = rep(pch, length.out = length(x))
  cex  = rep(cex, length.out = length(x))
  non_empty = !is.na(st_dimension(x))
  lapply(seq_along(x), function(i)
    if (non_empty[i])
      lines(x[[i]], lty = lty[i], lwd = lwd[i], col = col[i], pch = pch[i], cex = cex[i], type = type))
  invisible(NULL)
}

#' xy_beeline_len
#'
#' @param pt
#'
#' @return
#' @export
#'
#' @examples
st_beeline <- function(pt) {
  pt[c(1,length(pt))] %>% st_pt_line()
}

#' @export
st_pt_line <- function(pt)
  pt %>%
  do.call(rbind,.) %>%
  st_linestring() %>% #a single line
  st_sfc(crs=st_crs(pt))

#' st_azimuth
#'
#' description here ....
#'
#' @param x
#' @param y
#' @param ...
#'
#' @return
#' @export
#'
#' @examples
#' st_sfc(st_point(c(1,1)),
#'        st_point(c(3,3)),
#'        crs=get_crs('wgs84')) %>%
#' st_azimuth()
st_azimuth <- function(x,y,...) UseMethod("st_azimuth")
#fallback method
st_azimuth.default <- function(x,y,...) "st_azimuth is only implemented for sfc_POINT"

#' @export
st_azimuth.sfc_POINT <- function(x,y,bearingRhumb=F,...){
  if (missing(y))
    y = x
  else
    stopifnot(st_crs(x) == st_crs(y))
  x = st_geometry(x)
  y = st_geometry(y)
  if (!is.na(st_crs(x)))
    p = sf:::crs_parameters(st_crs(x)) #not exported

  if (isTRUE(st_is_longlat(x))) {
    if (!inherits(x, "sfc_POINT") || !inherits(y, "sfc_POINT")) #for multiple dispath
      stop("st_azimuth for longitude/latitude data only available for POINT geometries")
    if (!requireNamespace("geosphere", quietly = TRUE))
      stop("package geosphere required, please install it first")

    xp = do.call(rbind, x)[rep(seq_along(x), length(y)),]
    yp = do.call(rbind, y)[rep(seq_along(y), each = length(x)),]

    m = matrix(
      if (bearingRhumb)
        geosphere::bearingRhumb(xp,yp)
      else
        geosphere::bearing(xp, yp, as.numeric(p$SemiMajor), 1./p$InvFlattening),
      length(x), length(y))
    m <- set_units(m,degree)
    m
  } else {
    do.call(atan2, rev(as.list(st_point(c(1, 1)) -
                                 st_point(c(0,0)))
    )) * 180/pi

    # todo check geos and open an issue
    # where is it maintainted? opengeo
    # d = CPL_geos_dist(x, y)
    # if (! is.na(st_crs(x)))
    #   units(d) = p$ud_unit
    # d
  }
}

#' @export
st_segment <- function(x1,y1,x2,y2,crs=NA_crs_){
  li <- recycle_args(x1,y1,x2,y2)
  df <- as.data.frame(t(do.call(cbind,li)))
  linestrings <- lapply(df,
                        function(q)
                          sf::st_linestring(matrix(q,ncol=2,byrow=T),dim = "XY"))
  st_sfc(linestrings,crs=crs)
}

#' @export
st_sfc_point <- function(x,y,crs){
  #st_sf()
  mapply(x,y,
         FUN = function(x, y)
           st_point(c(x, y)),
         SIMPLIFY = F) %>%
    st_sfc(crs=crs)
}

#' st nearest dist
#'
#' @param sfc1
#' @param sfc2
#'
#' @return distance to nearest spatial object
#' @export
#'
#' @examples
st_nearest_dist <- function(sfc1,sfc2) {
  sf::st_distance(sfc1,sfc2) %>%
    apply(1,min)
}

#' read shapes
#' @description uses sf to read shapefiles of a folder
#' @param folder
#'
#' @return
#' @export
#'
#' @examples
st_read_shps <- function(folder){
  files <- list.files(folder,pattern = "*.shp$",full.names = T)
  shps <- lapply(files, sf::st_read, quiet=T)
  #regmatches(files,regexec(pat,files,perl=T))
  names(shps) <- gsub(pattern = ".+/(.+)\\.shp",
                      files,
                      replacement = "\\1")
  shps
}

#' st_drop_geometry
#'
#' @param x sf
#'
#' @return data.frame
#' @export
#'
#' @examples
st_drop_geometry <- function(x) {
  # drop=T also works
  # if(inherits(x,"sf"))
  #   ret <- x[,setdiff(names(x),attr(x,'sf_column')),drop=T]
  # else
  #   ret <- x
  #
  # class(ret) <- 'data.frame'
  # return(ret)

  if(inherits(x,"sf")) {
    x <- st_set_geometry(x, NULL)
    class(x) <- 'data.frame'
  }
  return(x)
}

#' Get Common CRS
#' @description gets the proj4 string of most useful crs by their short names.
#' google `epsg utm 39` for adding other crs.
#' @param crsName crs name
#'
#' @return the proj4 string is returned
#' @export
#'
#' @examples
#' get_crs("gmap")
st_get_crs <- function(crsName=c("wgs84","utm40","utm39","gearth","gmap","utm40_km")){

  formal.args <- formals(sys.function())

  #default is first
  crs <- match.arg(crsName)
  a <- c(
    "+proj=longlat +ellps=WGS84",
    "+init=epsg:32640",
    "+init=epsg:32639",
    "+init=epsg:4326",
    "+init=epsg:3857",
    "+proj=utm +zone=40 +datum=WGS84 +units=km +no_defs"
  )
  names(a) <- formal.args$crsName

  return(st_crs(a[crs]))
}