591 R/SpatialLinesNetwork.R
View file

Large diffs are not rendered by default.

137 R/aggregate_funs.R
View file
@@ -32,10 +32,10 @@
#'
#' @export
#' @examples
#' zones$quadrant = c(1, 2, 1, 4, 5, 6, 7, 1)
#' zones$quadrant <- c(1, 2, 1, 4, 5, 6, 7, 1)
#' aggzones <- rgeos::gUnaryUnion(zones, id = zones@data$quadrant)
#' aggzones <- sp::SpatialPolygonsDataFrame(aggzones, data.frame(region = c(1:6)), match.ID = FALSE)
#' sp::proj4string(aggzones) = sp::proj4string(zones)
#' sp::proj4string(aggzones) <- sp::proj4string(zones)
#' aggzones_sf <- sf::st_as_sf(aggzones)
#' aggzones_sf <- sf::st_set_crs(aggzones_sf, sf::st_crs(zones_sf))
#' od_agg <- od_aggregate(flow, zones_sf, aggzones_sf)
@@ -60,7 +60,6 @@ od_aggregate.sf <- function(flow, zones, aggzones,
FUN = sum,
prop_by_area = ifelse(identical(FUN, mean) == FALSE, TRUE, FALSE),
digits = getOption("digits")) {

flow_first_col <- colnames(flow)[1]
flow_second_col <- colnames(flow)[2]
zonesfirstcol <- colnames(zones)[1]
@@ -75,7 +74,7 @@ od_aggregate.sf <- function(flow, zones, aggzones,
}

zone_points <- sf::st_centroid(zones)
if(is.null(aggzone_points)) {
if (is.null(aggzone_points)) {
aggzone_points <- sf::st_centroid(aggzones)
}

@@ -84,8 +83,8 @@ od_aggregate.sf <- function(flow, zones, aggzones,
sf::st_set_geometry(NULL)

names(zones_agg)[1] <- flow_first_col
zones_agg$new_orig = zones_agg[, aggcols[1]]
zones_agg$new_dest = zones_agg[, aggcols[1]]
zones_agg$new_orig <- zones_agg[, aggcols[1]]
zones_agg$new_dest <- zones_agg[, aggcols[1]]

flow_new_orig <- flow %>%
dplyr::inner_join(y = zones_agg[c(flow_first_col, "new_orig")])
@@ -103,7 +102,6 @@ od_aggregate.sf <- function(flow, zones, aggzones,
flow_ag

# od2line(flow = flow_ag, zones = aggzones) # to export as sf

}
#' @export
od_aggregate.Spatial <- function(flow, zones, aggzones,
@@ -117,7 +115,7 @@ od_aggregate.Spatial <- function(flow, zones, aggzones,
aggzonesfirstcol <- colnames(aggzones@data)[1]

if (cols == FALSE) {
cols <- unlist(lapply(flow, is, 'numeric'))
cols <- unlist(lapply(flow, is, "numeric"))
cols <- names(cols[which(cols == TRUE)])
}
if (aggcols == FALSE) {
@@ -130,9 +128,11 @@ od_aggregate.Spatial <- function(flow, zones, aggzones,
if (sp::is.projected(zones) == TRUE & all.equal(zones@proj4string, aggzones@proj4string) == FALSE) {
aggzones <- sp::spTransform(aggzones, zones@proj4string)
} else {
projection <- paste0("+proj=aea +lat_1=90 +lat_2=-18.416667 ",
"+lat_0=0 +lon_0=10 +x_0=0 +y_0=0 +ellps=GRS80",
" +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
projection <- paste0(
"+proj=aea +lat_1=90 +lat_2=-18.416667 ",
"+lat_0=0 +lon_0=10 +x_0=0 +y_0=0 +ellps=GRS80",
" +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
)
zones <- sp::spTransform(zones, projection)
aggzones <- sp::spTransform(aggzones, projection)
}
@@ -143,50 +143,53 @@ od_aggregate.Spatial <- function(flow, zones, aggzones,

zoneintersect <- rgeos::gIntersection(zones, aggzones, byid = TRUE)
zoneintersect <- sp::SpatialPolygonsDataFrame(zoneintersect,
data=data.frame(
od_aggregate_charid=sapply(zoneintersect@polygons, function(x) x@ID),
row.names=sapply(zoneintersect@polygons, function(x) x@ID)
))
data = data.frame(
od_aggregate_charid = sapply(zoneintersect@polygons, function(x) x@ID),
row.names = sapply(zoneintersect@polygons, function(x) x@ID)
)
)
zoneintersect@data$od_aggregate_interarea <- rgeos::gArea(zoneintersect, byid = TRUE)
zoneintersect@data$od_aggregate_zone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[,1]
zoneintersect@data$od_aggregate_aggzone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[,2]
zoneintersect@data$od_aggregate_zone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[, 1]
zoneintersect@data$od_aggregate_aggzone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[, 2]

zoneintersect <- merge(zoneintersect, zones@data, by.x = 'od_aggregate_zone_charid', by.y = 'od_aggregate_charid')
zoneintersect@data$od_aggregate_proparea <- zoneintersect@data$od_aggregate_interarea/zoneintersect@data$stplanr_area
zoneintersect <- merge(zoneintersect, zones@data, by.x = "od_aggregate_zone_charid", by.y = "od_aggregate_charid")
zoneintersect@data$od_aggregate_proparea <- zoneintersect@data$od_aggregate_interarea / zoneintersect@data$stplanr_area

intersectdf <- merge(merge(
flow,
setNames(zoneintersect@data, paste0('o_', colnames(zoneintersect@data))),
by.x=colnames(flow)[1],
by.y=paste0('o_',zonesfirstcol)),
setNames(zoneintersect@data, paste0('d_', colnames(zoneintersect@data))),
by.x=colnames(flow)[2],
by.y=paste0('d_',zonesfirstcol)
setNames(zoneintersect@data, paste0("o_", colnames(zoneintersect@data))),
by.x = colnames(flow)[1],
by.y = paste0("o_", zonesfirstcol)
),
setNames(zoneintersect@data, paste0("d_", colnames(zoneintersect@data))),
by.x = colnames(flow)[2],
by.y = paste0("d_", zonesfirstcol)
)

if (prop_by_area == TRUE & is(zones, "SpatialPolygonsDataFrame") == TRUE) {
intersectdf <- intersectdf %>%
dplyr::mutate_at(
cols, dplyr::funs_('round(.*o_od_aggregate_proparea*d_od_aggregate_proparea)',args = list('digits'=digits))
cols, dplyr::funs_("round(.*o_od_aggregate_proparea*d_od_aggregate_proparea)", args = list("digits" = digits))
)
}

intersectdf <- intersectdf %>%
dplyr::group_by_('o_od_aggregate_aggzone_charid', 'd_od_aggregate_aggzone_charid') %>%
dplyr::select(dplyr::one_of(c('o_od_aggregate_aggzone_charid','d_od_aggregate_aggzone_charid',cols))) %>%
dplyr::summarise_at(cols,.funs = FUN) %>%
dplyr::left_join(setNames(aggzones@data[,c('od_aggregate_charid', aggcols)], c('od_aggregate_charid', paste0('o_',aggcols))),
by = c('o_od_aggregate_aggzone_charid' = 'od_aggregate_charid')) %>%
dplyr::left_join(setNames(aggzones@data[,c('od_aggregate_charid', aggcols)], c('od_aggregate_charid', paste0('d_',aggcols))),
by = c('d_od_aggregate_aggzone_charid' = 'od_aggregate_charid'))
intersectdf <- intersectdf[,c(
paste0('o_', c(aggzonesfirstcol, aggcols[which(aggcols != aggzonesfirstcol)])),
paste0('d_', c(aggzonesfirstcol, aggcols[which(aggcols != aggzonesfirstcol)])),
dplyr::group_by_("o_od_aggregate_aggzone_charid", "d_od_aggregate_aggzone_charid") %>%
dplyr::select(dplyr::one_of(c("o_od_aggregate_aggzone_charid", "d_od_aggregate_aggzone_charid", cols))) %>%
dplyr::summarise_at(cols, .funs = FUN) %>%
dplyr::left_join(setNames(aggzones@data[, c("od_aggregate_charid", aggcols)], c("od_aggregate_charid", paste0("o_", aggcols))),
by = c("o_od_aggregate_aggzone_charid" = "od_aggregate_charid")
) %>%
dplyr::left_join(setNames(aggzones@data[, c("od_aggregate_charid", aggcols)], c("od_aggregate_charid", paste0("d_", aggcols))),
by = c("d_od_aggregate_aggzone_charid" = "od_aggregate_charid")
)
intersectdf <- intersectdf[, c(
paste0("o_", c(aggzonesfirstcol, aggcols[which(aggcols != aggzonesfirstcol)])),
paste0("d_", c(aggzonesfirstcol, aggcols[which(aggcols != aggzonesfirstcol)])),
cols
)]

return(as.data.frame(intersectdf))

}

#' Aggregate SpatialPolygonsDataFrame to new geometry.
@@ -216,11 +219,11 @@ od_aggregate.Spatial <- function(flow, zones, aggzones,
#' @examples
#' \dontrun{
#' zones@data$region <- 1
#' zones@data[c(2, 5), c('region')] <- 2
#' zones@data[c(2, 5), c("region")] <- 2
#' aggzones <- sp::SpatialPolygonsDataFrame(rgeos::gUnaryUnion(
#' zones,
#' id = zones@data$region), data.frame(region=c(1, 2))
#' )
#' zones,
#' id = zones@data$region
#' ), data.frame(region = c(1, 2)))
#' zones@data$region <- NULL
#' zones@data$exdata <- 5
#' library(sp)
@@ -229,14 +232,13 @@ od_aggregate.Spatial <- function(flow, zones, aggzones,
sp_aggregate <- function(zones, aggzones, cols = FALSE,
FUN = sum,
prop_by_area = ifelse(identical(FUN, mean) == FALSE, TRUE, FALSE),
digits = getOption("digits")){

digits = getOption("digits")) {
zonesfirstcol <- colnames(zones@data)[1]
aggzonesfirstcol <- colnames(aggzones@data)[1]
aggcols <- colnames(aggzones@data)

if (cols == FALSE) {
cols <- unlist(lapply(zones@data, is, 'numeric'))
cols <- unlist(lapply(zones@data, is, "numeric"))
cols <- names(cols[which(cols == TRUE)])
cols <- cols[which(cols != zonesfirstcol)]
}
@@ -247,9 +249,11 @@ sp_aggregate <- function(zones, aggzones, cols = FALSE,
if (sp::is.projected(zones) == TRUE & all.equal(zones@proj4string, aggzones@proj4string) == FALSE) {
aggzones <- sp::spTransform(aggzones, zones@proj4string)
} else {
projection <- paste0("+proj=aea +lat_1=90 +lat_2=-18.416667 ",
"+lat_0=0 +lon_0=10 +x_0=0 +y_0=0 +ellps=GRS80",
" +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
projection <- paste0(
"+proj=aea +lat_1=90 +lat_2=-18.416667 ",
"+lat_0=0 +lon_0=10 +x_0=0 +y_0=0 +ellps=GRS80",
" +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
)
zones <- sp::spTransform(zones, projection)
aggzones <- sp::spTransform(aggzones, projection)
}
@@ -260,35 +264,39 @@ sp_aggregate <- function(zones, aggzones, cols = FALSE,

zoneintersect <- rgeos::gIntersection(zones, aggzones, byid = TRUE)
zoneintersect <- sp::SpatialPolygonsDataFrame(zoneintersect,
data=data.frame(
od_aggregate_charid=sapply(zoneintersect@polygons, function(x) x@ID),
row.names=sapply(zoneintersect@polygons, function(x) x@ID)
))
data = data.frame(
od_aggregate_charid = sapply(zoneintersect@polygons, function(x) x@ID),
row.names = sapply(zoneintersect@polygons, function(x) x@ID)
)
)
zoneintersect@data$od_aggregate_interarea <- rgeos::gArea(zoneintersect, byid = TRUE)
zoneintersect@data$od_aggregate_zone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[,1]
zoneintersect@data$od_aggregate_aggzone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[,2]
zoneintersect@data$od_aggregate_zone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[, 1]
zoneintersect@data$od_aggregate_aggzone_charid <- stringr::str_split(zoneintersect@data$od_aggregate_charid, " ", simplify = TRUE)[, 2]

zoneintersect <- merge(zoneintersect, zones@data, by.x = 'od_aggregate_zone_charid', by.y = 'od_aggregate_charid')
zoneintersect@data$od_aggregate_proparea <- zoneintersect@data$od_aggregate_interarea/zoneintersect@data$stplanr_area
zoneintersect <- merge(zoneintersect, zones@data, by.x = "od_aggregate_zone_charid", by.y = "od_aggregate_charid")
zoneintersect@data$od_aggregate_proparea <- zoneintersect@data$od_aggregate_interarea / zoneintersect@data$stplanr_area

intersectdf <- zoneintersect@data

if (prop_by_area == TRUE & is(zones, "SpatialPolygonsDataFrame") == TRUE) {
intersectdf <- intersectdf %>%
dplyr::mutate_at(
cols, dplyr::funs_('round(.*od_aggregate_proparea)',args = list('digits'=digits))
cols, dplyr::funs_("round(.*od_aggregate_proparea)", args = list("digits" = digits))
)
}

intersectdf <- intersectdf %>%
dplyr::group_by_('od_aggregate_aggzone_charid') %>%
dplyr::select(dplyr::one_of(c('od_aggregate_aggzone_charid',cols))) %>%
dplyr::summarise_at(cols,.funs = FUN) %>%
dplyr::left_join(setNames(aggzones@data[,c('od_aggregate_charid', aggcols)],c('od_aggregate_aggzone_charid', aggcols)),
by = 'od_aggregate_aggzone_charid')
intersectdf <- as.data.frame(intersectdf,
intersectdf$od_aggregate_aggzone_charid)
intersectdf <- intersectdf[,c(
dplyr::group_by_("od_aggregate_aggzone_charid") %>%
dplyr::select(dplyr::one_of(c("od_aggregate_aggzone_charid", cols))) %>%
dplyr::summarise_at(cols, .funs = FUN) %>%
dplyr::left_join(setNames(aggzones@data[, c("od_aggregate_charid", aggcols)], c("od_aggregate_aggzone_charid", aggcols)),
by = "od_aggregate_aggzone_charid"
)
intersectdf <- as.data.frame(
intersectdf,
intersectdf$od_aggregate_aggzone_charid
)
intersectdf <- intersectdf[, c(
c(aggzonesfirstcol, aggcols[which(aggcols != aggzonesfirstcol)]),
cols
)]
@@ -297,5 +305,4 @@ sp_aggregate <- function(zones, aggzones, cols = FALSE,
aggzones@data <- intersectdf

return(aggzones)

}
10 R/as_sf_fun.R
View file
@@ -5,23 +5,23 @@
#' @param ... Arguments passed to \code{FUN}
#' @aliases as_sp_fun
as_sf_fun <- function(input, FUN, ...) {
if(is(object = input, class2 = "sf")) {
if (is(object = input, class2 = "sf")) {
input <- as(object = input, Class = "Spatial")
}
res <- FUN(input)
if(is(object = res, class2 = "Spatial")) {
if (is(object = res, class2 = "Spatial")) {
res <- sf::st_as_sf(res)
}
return(res)
}

as_sp_fun <- function(input, FUN, ...) {
if(is(object = input, class2 = "Spatial")) {
if (is(object = input, class2 = "Spatial")) {
input <- sf::st_as_sf(input)
}
res <- FUN(input)
if(is(object = res, class2 = "sf")) {
if (is(object = res, class2 = "sf")) {
res <- as(res, "Spatial")
}
return(res)
}
}
445 R/catchmentArea.R
View file

Large diffs are not rendered by default.

26 R/crs-funs.R
View file
@@ -20,10 +20,12 @@
#' line_length <- rgeos::gLength(rf_projected, byid = TRUE)
#' plot(line_length, rf_projected$length)
#' cor(line_length, rf_projected$length)
crs_select_aeq <- function(shp){
crs_select_aeq <- function(shp) {
cent <- rgeos::gCentroid(shp)
aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
cent@coords[[2]], cent@coords[[1]])
aeqd <- sprintf(
"+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
cent@coords[[2]], cent@coords[[1]]
)
sp::CRS(aeqd)
}
#' Reproject lat/long spatial object so that they are in units of 1m
@@ -36,17 +38,17 @@ crs_select_aeq <- function(shp){
#' @export
#' @examples
#' data(routes_fast)
#' rf_aeq = reproject(routes_fast[1:3, ])
#' rf_osgb = reproject(routes_fast[1:3, ], 27700)
reproject <- function(shp, crs = crs_select_aeq(shp)){
if(is.na(raster::crs(shp))){
#' rf_aeq <- reproject(routes_fast[1:3, ])
#' rf_osgb <- reproject(routes_fast[1:3, ], 27700)
reproject <- function(shp, crs = crs_select_aeq(shp)) {
if (is.na(raster::crs(shp))) {
message("Assuming a geographical (lat/lon) CRS (EPSG:4326)")
raster::crs(shp) = sp::CRS("+init=epsg:4326")
raster::crs(shp) <- sp::CRS("+init=epsg:4326")
}
if (is.numeric(crs)) { # test if it's an epsg code
crs <- sp::CRS(paste0("+init=epsg:", crs))
}
if(is.numeric(crs)) # test if it's an epsg code
crs = sp::CRS(paste0("+init=epsg:", crs))
message(paste0("Transforming to CRS ", crs))
res = sp::spTransform(shp, crs)
res <- sp::spTransform(shp, crs)
res
}

22 R/cyclestreets.R
View file
@@ -12,22 +12,24 @@
#' coordinate to map.
#'
#' @export
#' @examples \dontrun{
#' nearest_cyclestreets(53, 0.02, pat = Sys.getenv("CYCLESTREET"))
#' nearest_cyclestreets(cents[1, ], pat = Sys.getenv("CYCLESTREET"))
#' nearest_cyclestreets(cents_sf[1, ], pat = Sys.getenv("CYCLESTREET"))
#' @examples
#' \dontrun{
#' nearest_cyclestreets(53, 0.02, pat = Sys.getenv("CYCLESTREET"))
#' nearest_cyclestreets(cents[1, ], pat = Sys.getenv("CYCLESTREET"))
#' nearest_cyclestreets(cents_sf[1, ], pat = Sys.getenv("CYCLESTREET"))
#' }
nearest_cyclestreets <- function(shp = NULL, lat, lng, pat = api_pat("cyclestreet")) {
UseMethod("nearest_cyclestreets", object = shp)
}
#' @export
nearest_cyclestreets.NULL <- function(shp = NULL, lat, lng, pat = api_pat("cyclestreet")) {
url = paste0("https://api.cyclestreets.net/v2/nearestpoint?lonlat=", lng, ",", lat, "&key=", pat)
obj = jsonlite::fromJSON(url)
coords = obj$features$geometry$coordinates[[1]]
sp::SpatialPointsDataFrame(coords = matrix(coords, ncol = 2),
data = data.frame(orig_lat = lat, orig_lng = lng))

url <- paste0("https://api.cyclestreets.net/v2/nearestpoint?lonlat=", lng, ",", lat, "&key=", pat)
obj <- jsonlite::fromJSON(url)
coords <- obj$features$geometry$coordinates[[1]]
sp::SpatialPointsDataFrame(
coords = matrix(coords, ncol = 2),
data = data.frame(orig_lat = lat, orig_lng = lng)
)
}
#' @export
nearest_cyclestreets.Spatial <- function(shp, lat = shp@coords[1, 2], lng = shp@coords[1, 1], pat = api_pat("cyclestreet")) {
65 R/data.r
View file
@@ -24,17 +24,17 @@
#' home <- sp::spTransform(x = home, CRSobj = crsuk)
#' buf <- rgeos::gBuffer(home, width = 2000)
#' # Check it saved the points OK
#' cents <- cents[buf,]
#' cents <- cents[buf, ]
#' plot(buf)
#' points(cents)
#' cents <- sp::spTransform(x = cents, CRSobj = crs)
#' cents$geo_code <- as.character(cents$geo_code)
#' library(devtools)
#' # use_data(cents, overwrite = TRUE)
#' cents_sf = sf::st_as_sf(cents)
#' cents_sf <- sf::st_as_sf(cents)
#' devtools::use_data(cents_sf)
#' }
#'
#'
#' @docType data
#' @keywords datasets
#' @name cents
@@ -76,21 +76,21 @@ NULL
#' library(devtools)
#' flow$id <- paste(flow$Area.of.residence, flow$Area.of.workplace)
#' use_data(flow, overwrite = TRUE)
#'
#'
#' # Convert flows to spatial lines dataset
#' flowlines <- od2line(flow = flow, zones = cents)
#' # use_data(flowlines, overwrite = TRUE)
#'
#'
#' # Convert flows to routes
#' routes_fast <- line2route(l = flowlines, plan = "fastest")
#' routes_slow <- line2route(l = flowlines, plan = "quietest")
#'
#'
#' use_data(routes_fast)
#' use_data(routes_slow)
#' routes_fast_sf <- sf::st_as_sf(routes_fast)
#' routes_slow_sf <- sf::st_as_sf(routes_slow)
#' }
#'
#'
#' @docType data
#' @keywords datasets
#' @name flow
@@ -105,12 +105,12 @@ NULL
#' @examples
#' \dontrun{
#' # This is how the dataset was constructed
#' flow_dests = flow
#' flow_dests$Area.of.workplace = sample(x = destinations$WZ11CD, size = nrow(flow))
#' flow_dests = dplyr::rename(flow_dests, WZ11CD = Area.of.workplace)
#' flow_dests <- flow
#' flow_dests$Area.of.workplace <- sample(x = destinations$WZ11CD, size = nrow(flow))
#' flow_dests <- dplyr::rename(flow_dests, WZ11CD = Area.of.workplace)
#' devtools::use_data(flow_dests)
#' }
#'
#'
#' @docType data
#' @keywords datasets
#' @name flow_dests
@@ -127,22 +127,24 @@ NULL
#' \dontrun{
#' # This is how the dataset was constructed - see
#' # http://cowz.geodata.soton.ac.uk/download/
#' download.file("http://cowz.geodata.soton.ac.uk/download/files/COWZ_EW_2011_BFC.zip",
#' "COWZ_EW_2011_BFC.zip")
#' download.file(
#' "http://cowz.geodata.soton.ac.uk/download/files/COWZ_EW_2011_BFC.zip",
#' "COWZ_EW_2011_BFC.zip"
#' )
#' unzip("COWZ_EW_2011_BFC.zip")
#' wz = raster::shapefile("COWZ_EW_2011_BFC.shp")
#' to_remove = list.files(pattern = "COWZ", full.names = TRUE, recursive = TRUE)
#' wz <- raster::shapefile("COWZ_EW_2011_BFC.shp")
#' to_remove <- list.files(pattern = "COWZ", full.names = TRUE, recursive = TRUE)
#' file.remove(to_remove)
#' proj4string(wz)
#' wz = sp::spTransform(wz, proj4string(zones))
#' destination_zones = wz[zones,]
#' wz <- sp::spTransform(wz, proj4string(zones))
#' destination_zones <- wz[zones, ]
#' plot(destination_zones)
#' devtools::use_data(destination_zones)
#' head(destination_zones@data)
#' destinations = rgeos::gCentroid(destinations, byid = TRUE)
#' destinations = sp::SpatialPointsDataFrame(destinations, destination_zones@data)
#' destinations <- rgeos::gCentroid(destinations, byid = TRUE)
#' destinations <- sp::SpatialPointsDataFrame(destinations, destination_zones@data)
#' devtools::use_data(destinations, overwrite = TRUE)
#' destinations_sf = sf::st_as_sf(destinations)
#' destinations_sf <- sf::st_as_sf(destinations)
#' devtools::use_data(destinations_sf)
#' }
#' @docType data
@@ -212,10 +214,10 @@ NULL
#' \dontrun{
#' zones <- rgdal::readOGR(dsn = "/home/robin/npct/pct-bigdata/msoas.geojson", layer = "OGRGeoJSON")
#' proj4string(zones) <- proj4string(cents)
#' zones <- zones[cents,]
#' zones <- zones[cents, ]
#' plot(zones)
#' points(cents)
#' zones_sf = sf::st_as_sf(zones)
#' zones_sf <- sf::st_as_sf(zones)
#' }
#' @docType data
#' @keywords datasets
@@ -236,9 +238,10 @@ NULL
#' @aliases route_network_sf
#' @usage data(route_network)
#' @format A spatial lines dataset 80 rows and 1 column
#' @examples \dontrun{
#' @examples
#' \dontrun{
#' # Generate route network
#' route_network = overline(routes_fast, "All", fun = sum)
#' route_network <- overline(routes_fast, "All", fun = sum)
#' route_network_sf <- sf::st_as_sf(route_network)
#' }
NULL
@@ -254,7 +257,8 @@ NULL
#' @name ca_local
#' @usage data(ca_local)
#' @format A SpatialPointsDataFrame with 11 rows and 2 columns
#' @examples \dontrun{
#' @examples
#' \dontrun{
#' # Generate data
#' ac <- read_stats19_ac()
#' ca <- read_stats19_ca()
@@ -268,7 +272,7 @@ NULL
#' data("route_network")
#' proj4string(ca_sp) <- proj4string(route_network)
#' bb <- bb2poly(route_network)
#' ca_local = ca_sp[bb,]
#' ca_local <- ca_sp[bb, ]
#' }
NULL

@@ -281,12 +285,13 @@ NULL
#' @usage data(l_poly)
#' @format A SpatialPolygon
#'
#' @examples \dontrun{
#' l = routes_fast[13,]
#' l_poly = buff_geo(l, 8)
#' @examples
#' \dontrun{
#' l <- routes_fast[13, ]
#' l_poly <- buff_geo(l, 8)
#' plot(l_poly)
#' plot(routes_fast, add = TRUE)
#' # allocate road width to relevant line
#' devtools::use_data(l_poly)
#' }
NULL
NULL
104 R/geo-functions.R
View file
@@ -6,8 +6,8 @@
#'
#' @inheritParams gclip
#' @param filename File name of the output geojson
writeGeoJSON <- function(shp, filename){
name <- nm <-deparse(substitute(shp))
writeGeoJSON <- function(shp, filename) {
name <- nm <- deparse(substitute(shp))
rgdal::writeOGR(obj = shp, layer = name, dsn = filename, driver = "GeoJSON")
newname <- paste0(filename, ".geojson")
file.rename(filename, newname)
@@ -44,37 +44,37 @@ writeGeoJSON <- function(shp, filename){
#' @export
#' @examples
#' \dontrun{
#' shp <- routes_fast[2,]
#' shp <- routes_fast[2, ]
#' plot(shp)
#' rfs10 <- mapshape(shp)
#' rfs5 <- mapshape(shp, percent = 5)
#' rfs1 <- mapshape(shp, percent = 1)
#' plot(rfs10, add = TRUE, col ="red")
#' plot(rfs5, add = TRUE, col ="blue")
#' plot(rfs10, add = TRUE, col = "red")
#' plot(rfs5, add = TRUE, col = "blue")
#' plot(rfs1, add = TRUE, col = "grey")
#' # snap the lines to the nearest interval
#' rfs_int <- mapshape(shp, ms_options = "snap-interval=0.001")
#' plot(shp)
#' plot(rfs_int, add = TRUE)
#' mapshape(routes_fast_sf[2,])
#' mapshape(routes_fast_sf[2, ])
#' }
mapshape <- function(shp, percent = 10, ms_options = "", dsn = "mapshape", silent = FALSE){
shp_filename = paste0(dsn, ".shp")
new_filename = paste0(dsn, "-ms.shp")
if(!mapshape_available()) stop("mapshaper not available on this system")
mapshape <- function(shp, percent = 10, ms_options = "", dsn = "mapshape", silent = FALSE) {
shp_filename <- paste0(dsn, ".shp")
new_filename <- paste0(dsn, "-ms.shp")
if (!mapshape_available()) stop("mapshaper not available on this system")
is_sp <- is(shp, "Spatial")
if(is_sp) {
if (is_sp) {
shp <- sf::st_as_sf(shp)
}
sf::write_sf(shp, shp_filename, delete_layer = TRUE)
cmd <- paste0("mapshaper ", ms_options, " ", shp_filename, " -simplify ", percent, "% -o ", new_filename)
if(!silent) print(paste0("Running the following command from the system: ", cmd))
if (!silent) print(paste0("Running the following command from the system: ", cmd))
system(cmd, ignore.stderr = TRUE)
new_shp <- sf::st_read(paste0(dsn, "-ms.shp"))
sf::st_crs(new_shp) <- sf::st_crs(shp)
to_remove <- list.files(pattern = dsn)
file.remove(to_remove)
if(is_sp) {
if (is_sp) {
new_shp <- as(new_shp, "Spatial")
}
new_shp
@@ -114,22 +114,22 @@ mapshape_available <- function() {
#' plot(clipped, add = TRUE)
#' clipped$avslope # gclip also returns the data attribute
#' points(clipped)
#' points(cents[cb,], col = "red") # note difference
#' points(cents[cb, ], col = "red") # note difference
#' gclip(cents_sf, cb)
gclip <- function(shp, bb) {
UseMethod("gclip")
}
#' @export
gclip.Spatial <- function(shp, bb) {
if(class(bb) == "matrix"){
if (class(bb) == "matrix") {
b_poly <- as(raster::extent(as.vector(t(bb))), "SpatialPolygons")
}
else{
else {
b_poly <- as(raster::extent(bb), "SpatialPolygons")
}
clipped <- rgeos::gIntersection(shp, b_poly, byid = TRUE, id = row.names(shp))
if(grepl("DataFrame", class(shp))){
if(grepl("SpatialLines", class(shp)) & grepl("SpatialCollections",class(clipped))) {
if (grepl("DataFrame", class(shp))) {
if (grepl("SpatialLines", class(shp)) & grepl("SpatialCollections", class(clipped))) {
geodata <- data.frame(gclip_id = row.names(clipped@lineobj))
}
else {
@@ -138,15 +138,15 @@ gclip.Spatial <- function(shp, bb) {
joindata <- cbind(gclip_id = row.names(shp), shp@data)
geodata <- dplyr::left_join(geodata, joindata)
row.names(geodata) <- geodata$gclip_id
#if the data are SpatialPolygonsDataFrame (based on https://stat.ethz.ch/pipermail/r-sig-geo/2008-January/003052.html)
if(grepl("SpatialPolygons", class(shp))){
#then rebuild SpatialPolygonsDataFrame selecting relevant rows by row.names (row ID values)
clipped <- sp::SpatialPolygonsDataFrame(clipped, as(shp[row.names(clipped),], "data.frame"))
} else if(grepl("SpatialLines", class(shp)) & grepl("SpatialCollections",class(clipped))) {
# if the data are SpatialPolygonsDataFrame (based on https://stat.ethz.ch/pipermail/r-sig-geo/2008-January/003052.html)
if (grepl("SpatialPolygons", class(shp))) {
# then rebuild SpatialPolygonsDataFrame selecting relevant rows by row.names (row ID values)
clipped <- sp::SpatialPolygonsDataFrame(clipped, as(shp[row.names(clipped), ], "data.frame"))
} else if (grepl("SpatialLines", class(shp)) & grepl("SpatialCollections", class(clipped))) {
clipped <- sp::SpatialLinesDataFrame(clipped@lineobj, geodata)
} else if(grepl("SpatialLines", class(shp))) {
} else if (grepl("SpatialLines", class(shp))) {
clipped <- sp::SpatialLinesDataFrame(clipped, geodata)
} else { #assumes the data is a SpatialPointsDataFrame
} else { # assumes the data is a SpatialPointsDataFrame
clipped <- sp::SpatialPointsDataFrame(clipped, geodata)
}
}
@@ -155,9 +155,9 @@ gclip.Spatial <- function(shp, bb) {
}
#' @export
gclip.sf <- function(shp, bb) {
shp <- as(shp, "Spatial")
shp <- gclip.Spatial(shp, as(bb, "Spatial"))
sf::st_as_sf(shp)
shp <- as(shp, "Spatial")
shp <- gclip.Spatial(shp, as(bb, "Spatial"))
sf::st_as_sf(shp)
}
#' Scale a bounding box
#'
@@ -173,12 +173,12 @@ gclip.sf <- function(shp, bb) {
#' bb1 <- bbox_scale(bb, scale_factor = 1.05)
#' bb2 <- bbox_scale(bb, scale_factor = c(2, 1.05))
#' bb3 <- bbox_scale(bb, 0.1)
#' plot(x = bb2[1,], y = bb2[2,])
#' points(bb1[1,], bb1[2,])
#' points(bb3[1,], bb3[2,])
#' points(bb[1,], bb[2,], col = "red")
bbox_scale <- function(bb, scale_factor){
if(length(scale_factor == 1)) scale_factor <- rep(scale_factor, 2)
#' plot(x = bb2[1, ], y = bb2[2, ])
#' points(bb1[1, ], bb1[2, ])
#' points(bb3[1, ], bb3[2, ])
#' points(bb[1, ], bb[2, ], col = "red")
bbox_scale <- function(bb, scale_factor) {
if (length(scale_factor == 1)) scale_factor <- rep(scale_factor, 2)
b <- (bb - rowMeans(bb)) * scale_factor + rowMeans(bb)
b
}
@@ -220,13 +220,13 @@ geo_bb.Spatial <- function(shp, scale_factor = 1, distance = 0, output = c("poly
bb <- bbox_scale(bb = bb, scale_factor = scale_factor)
bb <- bb2poly(bb = bb, distance = distance)
sp::proj4string(bb) <- sp::proj4string(shp)
if(output == "polygon") {
if (output == "polygon") {
return(bb)
} else if(output == "points") {
} else if (output == "points") {
bb_point <- sp::SpatialPoints(raster::geom(bb)[1:4, c(5, 6)])
sp::proj4string(bb_point) <- sp::proj4string(shp)
return(bb_point)
} else if(output == "bb") {
} else if (output == "bb") {
return(geo_bb_matrix(bb))
}
}
@@ -239,14 +239,14 @@ geo_bb.sf <- function(shp, scale_factor = 1, distance = 0, output = c("polygon",
bb_sp <- bb2poly(bb = bb, distance = distance)
bb <- sf::st_as_sf(bb_sp)
sf::st_crs(bb) <- sf::st_crs(shp)
if(output == "polygon") {
if (output == "polygon") {
return(bb)
} else if(output == "points") {
} else if (output == "points") {
bb_point <- sp::SpatialPoints(raster::geom(bb_sp)[1:4, c(5, 6)])
bb_point <- sf::st_as_sf(bb_point)
sf::st_crs(bb_point) = sf::st_crs(shp)
sf::st_crs(bb_point) <- sf::st_crs(shp)
return(bb_point)
} else if(output == "bb") {
} else if (output == "bb") {
return(geo_bb_matrix(bb))
}
}
@@ -259,47 +259,47 @@ geo_bb.bbox <- function(shp, scale_factor = 1, distance = 0, output = c("polygon
bb_sp <- bb2poly(bb = bb, distance = distance)
bb <- sf::st_as_sf(bb_sp)
sf::st_crs(bb) <- sf::st_crs(shp)
if(output == "polygon") {
if (output == "polygon") {
return(bb)
} else if(output == "points") {
} else if (output == "points") {
bb_point <- sp::SpatialPoints(raster::geom(bb_sp)[1:4, c(5, 6)])
bb_point <- sf::st_as_sf(bb_point)
sf::st_crs(bb_point) = sf::st_crs(shp)
sf::st_crs(bb_point) <- sf::st_crs(shp)
return(bb_point)
} else if(output == "bb") {
} else if (output == "bb") {
return(geo_bb_matrix(bb))
}
}

#' @export
geo_bb.matrix <- function(shp, scale_factor = 1, distance = 0, output = c("polygon", "points", "bb")) {
output <- match.arg(output)
if(nrow(shp) != 2) {
if (nrow(shp) != 2) {
bb <- geo_bb_matrix(shp)
} else {
bb <- shp
}
bb <- bbox_scale(bb = bb, scale_factor = scale_factor)
bb <- bb2poly(bb = bb, distance = distance)
if(output == "polygon") {
if (output == "polygon") {
return(bb)
} else if(output == "points") {
} else if (output == "points") {
bb_point <- sp::SpatialPoints(raster::geom(bb)[1:4, c(5, 6)])
return(bb_point)
} else if(output == "bb") {
} else if (output == "bb") {
return(geo_bb_matrix(bb))
}
}

#' @export
bb2poly <- function(bb, distance = 0){
if(class(bb) == "matrix"){
bb2poly <- function(bb, distance = 0) {
if (class(bb) == "matrix") {
b_poly <- as(raster::extent(as.vector(t(bb))), "SpatialPolygons")
} else {
b_poly <- as(raster::extent(bb), "SpatialPolygons")
proj4string(b_poly) <- proj4string(bb)
}
if(distance > 0) {
if (distance > 0) {
b_poly_buff <- geo_buffer(shp = b_poly, width = distance)
b_poly <- bb2poly(b_poly_buff)
}
64 R/geo_code.R
View file
@@ -19,51 +19,49 @@
#' }
geo_code <- function(address,
service = "nominatim",
base_url = "https://maps.google.com/maps/api/geocode/json",
return_all = FALSE,
pat = NULL
) {

if(service == "nominatim") {
if(base_url == "https://maps.google.com/maps/api/geocode/json") {
base_url <- "https://nominatim.openstreetmap.org"
base_url = "https://maps.google.com/maps/api/geocode/json",
return_all = FALSE,
pat = NULL) {
if (service == "nominatim") {
if (base_url == "https://maps.google.com/maps/api/geocode/json") {
base_url <- "https://nominatim.openstreetmap.org"
}
place_name <- address
query <- list(q = place_name, format = "json")
if (!return_all) {
query <- c(query, limit = 1)
}
place_name <- address
query <- list (q = place_name, format = "json")
if(!return_all) {
query <- c(query, limit = 1)
}

q_url <- httr::modify_url(base_url, query = query)
res <- httr::GET (q_url)
txt <- httr::content(res, as = "text", encoding = "UTF-8",
type = "application/xml")
obj <- jsonlite::fromJSON(txt)
res_df <- data.frame(lon = obj$lon, lat = obj$lat, name = obj$display_name)
lon_lat <- as.numeric(c(lon = obj$lon[1], lat = obj$lat[1]))
q_url <- httr::modify_url(base_url, query = query)
res <- httr::GET(q_url)
txt <- httr::content(res,
as = "text", encoding = "UTF-8",
type = "application/xml"
)
obj <- jsonlite::fromJSON(txt)
res_df <- data.frame(lon = obj$lon, lat = obj$lat, name = obj$display_name)
lon_lat <- as.numeric(c(lon = obj$lon[1], lat = obj$lat[1]))
} else {

query <- list(address = address, sensor = "false")
if(!is.null(pat)) {
if (!is.null(pat)) {
query <- c(query, key = pat)
}
u <- httr::modify_url(base_url, query = query)
res <- jsonlite::fromJSON(u)
if(res$status == "OVER_QUERY_LIMIT") {
stop(res$error_message)
}
u <- httr::modify_url(base_url, query = query)
res <- jsonlite::fromJSON(u)
if (res$status == "OVER_QUERY_LIMIT") {
stop(res$error_message)
}

res_df <- jsonlite::flatten(res$results)
lon_lat <- c(
lon = res_df$geometry.location.lng,
lat = res_df$geometry.location.lat
res_df <- jsonlite::flatten(res$results)
lon_lat <- c(
lon = res_df$geometry.location.lng,
lat = res_df$geometry.location.lat
)
}

if(return_all) {
if (return_all) {
return(res_df)
} else {
return(lon_lat)
}

}
57 R/geo_projected.R
View file
@@ -23,17 +23,19 @@ geo_select_aeq <- function(shp) {
UseMethod("geo_select_aeq")
}
#' @export
geo_select_aeq.Spatial <- function(shp){
geo_select_aeq.Spatial <- function(shp) {
crs_select_aeq(shp)
}
#' @export
geo_select_aeq.sf <- function(shp){
geo_select_aeq.sf <- function(shp) {
cent <- sf::st_geometry(shp)
coords <- sf::st_coordinates(shp)
coords_mat <- matrix(coords[,1:2], ncol = 2)
coords_mat <- matrix(coords[, 1:2], ncol = 2)
midpoint <- apply(coords_mat, 2, mean)
aeqd <- sprintf("+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
midpoint[2], midpoint[1])
aeqd <- sprintf(
"+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0",
midpoint[2], midpoint[1]
)
sf::st_crs(aeqd)
}

@@ -50,54 +52,57 @@ geo_select_aeq.sf <- function(shp){
#' @aliases gprojected
#' @export
#' @examples
#' shp = routes_fast_sf[2:4,]
#' shp <- routes_fast_sf[2:4, ]
#' plot(geo_projected(shp, sf::st_buffer, dist = 100)$geometry)
#' shp = routes_fast[2:4,]
#' shp <- routes_fast[2:4, ]
#' geo_projected(shp, fun = rgeos::gBuffer, width = 100, byid = TRUE)
#' rlength = geo_projected(routes_fast, fun = rgeos::gLength, byid = TRUE)
#' rlength <- geo_projected(routes_fast, fun = rgeos::gLength, byid = TRUE)
#' plot(routes_fast$length, rlength)
geo_projected <- function(shp, fun, crs, silent, ...) {
UseMethod(generic = "geo_projected")
}
#' @export
geo_projected.sf <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE, ...){
geo_projected.sf <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE, ...) {
# assume it's not projected (i.e. lat/lon) if there is no CRS
if(is.na(sf::st_crs(shp))) {
if (is.na(sf::st_crs(shp))) {
sf::st_crs(shp) <- 4326
}
crs_orig <- sf::st_crs(shp)
shp_projected <- sf::st_transform(shp, crs)
if(!silent) {
if (!silent) {
message(paste0("Running function on a temporary projection: ", crs$proj4string))
}
res <- fun(shp_projected, ...)
if(grepl("sf", x = class(res)[1]))
if (grepl("sf", x = class(res)[1])) {
res <- sf::st_transform(res, crs_orig)
}
res
}
#' @export
geo_projected.Spatial <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE, ...){
geo_projected.Spatial <- function(shp, fun, crs = geo_select_aeq(shp), silent = TRUE, ...) {
# assume it's not projected (i.e. lat/lon) if there is no CRS
if(!is.na(sp::is.projected(shp))){
if(sp::is.projected(shp)){
if (!is.na(sp::is.projected(shp))) {
if (sp::is.projected(shp)) {
res <- fun(shp, ...)
} else {
shp_projected <- reproject(shp, crs = crs)
if(!silent) {
if (!silent) {
message(paste0("Running function on a temporary projection: ", crs))
}
res <- fun(shp_projected, ...)
if(is(res, "Spatial"))
if (is(res, "Spatial")) {
res <- sp::spTransform(res, sp::CRS("+init=epsg:4326"))
}
}
} else {
shp_projected <- reproject(shp, crs = crs)
if(!silent) {
if (!silent) {
message(paste0("Running function on a temporary projection: ", crs$proj4string))
}
res <- fun(shp_projected, ...)
if(is(res, "Spatial"))
if (is(res, "Spatial")) {
res <- sp::spTransform(res, sp::CRS("+init=epsg:4326"))
}
}
res
}
@@ -114,11 +119,11 @@ gprojected <- geo_projected.Spatial
#' @param ... Arguments passed to the buffer (see \code{?rgeos::gBuffer} or \code{?sf::st_buffer} for details)
#' @seealso buff_geo
#' @examples
#' buff_sp = geo_buffer(routes_fast, width = 100)
#' buff_sp <- geo_buffer(routes_fast, width = 100)
#' class(buff_sp)
#' plot(buff_sp, col = "red")
#' routes_fast_sf = sf::st_as_sf(routes_fast)
#' buff_sf = geo_buffer(routes_fast_sf, dist = 50)
#' routes_fast_sf <- sf::st_as_sf(routes_fast)
#' buff_sf <- geo_buffer(routes_fast_sf, dist = 50)
#' plot(buff_sf$geometry, add = TRUE)
#' @export
geo_buffer <- function(shp, dist = NULL, width = NULL, ...) {
@@ -150,9 +155,11 @@ geo_length <- function(shp) {

#' @export
geo_length.sf <- function(shp) {
if(sf::st_is_longlat(shp))
l <- lwgeom::st_geod_length(shp) else
l <- sf::st_length(shp)
if (sf::st_is_longlat(shp)) {
l <- lwgeom::st_geod_length(shp)
} else {
l <- sf::st_length(shp)
}
as.numeric(l)
}

192 R/google-functions.R
View file
@@ -11,16 +11,19 @@
#' coordinate to map.
#' @param google_api String value containing the Google API key to use.
#' @export
#' @examples \dontrun{
#' nearest_google(lat = 50.333, lng = 3.222, google_api = "api_key_here")
#' @examples
#' \dontrun{
#' nearest_google(lat = 50.333, lng = 3.222, google_api = "api_key_here")
#' }
nearest_google <- function(lat, lng, google_api){
base_url = "https://roads.googleapis.com/v1/snapToRoads"
url = paste0(base_url, "?path=", lat, ",", lng, "&key=", google_api)
obj = jsonlite::fromJSON(url)
coords = c(obj$snappedPoints$location$longitude, obj$snappedPoints$location$latitude)
sp::SpatialPointsDataFrame(coords = matrix(coords, ncol = 2),
data = data.frame(orig_lat = lat, orig_lng = lng))
nearest_google <- function(lat, lng, google_api) {
base_url <- "https://roads.googleapis.com/v1/snapToRoads"
url <- paste0(base_url, "?path=", lat, ",", lng, "&key=", google_api)
obj <- jsonlite::fromJSON(url)
coords <- c(obj$snappedPoints$location$longitude, obj$snappedPoints$location$latitude)
sp::SpatialPointsDataFrame(
coords = matrix(coords, ncol = 2),
data = data.frame(orig_lat = lat, orig_lng = lng)
)
}
#' Return travel network distances and time using the Google Maps API
#'
@@ -43,110 +46,131 @@ nearest_google <- function(lat, lng, google_api){
#' simultaneous queries, and so will, for example, only returns values for up to
#' 10 origins times 10 destinations.
#'
#' @examples \dontrun{
#' # Distances from one origin to one destination
#' from = c(-46.3, -23.4)
#' to = c(-46.4, -23.4)
#' dist_google(from = from, to = to, mode = "walking") # not supported on last test
#' dist_google(from = from, to = to, mode = "driving")
#' dist_google(from = c(0, 52), to = c(0, 53))
#' data("cents")
#' # Distances from between all origins and destinations
#' dists_cycle = dist_google(from = cents, to = cents)
#' dists_drive = dist_google(cents, cents, mode = "driving")
#' dists_trans = dist_google(cents, cents, mode = "transit")
#' dists_trans_am = dist_google(cents, cents, mode = "transit",
#' @examples
#' \dontrun{
#' # Distances from one origin to one destination
#' from <- c(-46.3, -23.4)
#' to <- c(-46.4, -23.4)
#' dist_google(from = from, to = to, mode = "walking") # not supported on last test
#' dist_google(from = from, to = to, mode = "driving")
#' dist_google(from = c(0, 52), to = c(0, 53))
#' data("cents")
#' # Distances from between all origins and destinations
#' dists_cycle <- dist_google(from = cents, to = cents)
#' dists_drive <- dist_google(cents, cents, mode = "driving")
#' dists_trans <- dist_google(cents, cents, mode = "transit")
#' dists_trans_am <- dist_google(cents, cents,
#' mode = "transit",
#' arrival_time = strptime("2016-05-27 09:00:00",
#' format = "%Y-%m-%d %H:%M:%S", tz = "BST"))
#' # Find out how much longer (or shorter) cycling takes than walking
#' summary(dists_cycle$duration / dists_trans$duration)
#' # Difference between travelling now and for 9am arrival
#' summary(dists_trans_am$duration / dists_trans$duration)
#' odf = points2odf(cents)
#' odf = cbind(odf, dists)
#' head(odf)
#' flow = points2flow(cents)
#' # show the results for duration (thicker line = shorter)
#' plot(flow, lwd = mean(odf$duration) / odf$duration)
#' dist_google(c("Hereford"), c("Weobley", "Leominster", "Kington"))
#' dist_google(c("Hereford"), c("Weobley", "Leominster", "Kington"),
#' format = "%Y-%m-%d %H:%M:%S", tz = "BST"
#' )
#' )
#' # Find out how much longer (or shorter) cycling takes than walking
#' summary(dists_cycle$duration / dists_trans$duration)
#' # Difference between travelling now and for 9am arrival
#' summary(dists_trans_am$duration / dists_trans$duration)
#' odf <- points2odf(cents)
#' odf <- cbind(odf, dists)
#' head(odf)
#' flow <- points2flow(cents)
#' # show the results for duration (thicker line = shorter)
#' plot(flow, lwd = mean(odf$duration) / odf$duration)
#' dist_google(c("Hereford"), c("Weobley", "Leominster", "Kington"))
#' dist_google(c("Hereford"), c("Weobley", "Leominster", "Kington"),
#' mode = "transit", arrival_time = strptime("2016-05-27 17:30:00",
#' format = "%Y-%m-%d %H:%M:%S", tz = "BST"))
#' format = "%Y-%m-%d %H:%M:%S", tz = "BST"
#' )
#' )
#' }
dist_google <- function(from, to, google_api = Sys.getenv("GOOGLEDIST"),
g_units = 'metric',
g_units = "metric",
mode = c("bicycling", "walking", "driving", "transit"),
arrival_time = ""){
mode = match.arg(mode)
arrival_time = "") {
mode <- match.arg(mode)
base_url <- "https://maps.googleapis.com/maps/api/distancematrix/json?units="
# Convert sp object to lat/lon vector
if(class(from) == "SpatialPoints" | class(from) == "SpatialPointsDataFrame" )
if (class(from) == "SpatialPoints" | class(from) == "SpatialPointsDataFrame") {
from <- coordinates(from)
if(class(to) == "SpatialPoints" | class(to) == "SpatialPointsDataFrame" )
}
if (class(to) == "SpatialPoints" | class(to) == "SpatialPointsDataFrame") {
to <- coordinates(to)
}
if (google_api == "") {
google_api_param <- ""
} else {
google_api_param <- "&key="
}
if(is(from, "matrix") | is(from, "data.frame"))
from = paste(from[,2], from[,1], sep = ",")
if(is(from, "numeric"))
from = paste(from[2], from[1], sep = ",")
if(is(to, "matrix") | is(to, "data.frame"))
to = paste(to[,2], to[,1], sep = ",")
if(is(to, "numeric"))
to = paste(to[2], to[1], sep = ",")
from = paste0(from, collapse = "|")
to = paste0(to, collapse = "|")
url_travel <- paste0(base_url, g_units, "&origins=", from,
"&destinations=", to, "&mode=", mode)
if(class(arrival_time)[1] == "POSIXlt"){
if (is(from, "matrix") | is(from, "data.frame")) {
from <- paste(from[, 2], from[, 1], sep = ",")
}
if (is(from, "numeric")) {
from <- paste(from[2], from[1], sep = ",")
}
if (is(to, "matrix") | is(to, "data.frame")) {
to <- paste(to[, 2], to[, 1], sep = ",")
}
if (is(to, "numeric")) {
to <- paste(to[2], to[1], sep = ",")
}
from <- paste0(from, collapse = "|")
to <- paste0(to, collapse = "|")
url_travel <- paste0(
base_url, g_units, "&origins=", from,
"&destinations=", to, "&mode=", mode
)
if (class(arrival_time)[1] == "POSIXlt") {
arrival_time <- as.numeric(arrival_time)
url_travel <- paste0(url_travel, "&arrival_time=", arrival_time)
}
url = paste0(url_travel,
google_api_param,
google_api)
url = utils::URLencode(url, repeated = FALSE, reserved = FALSE)
url <- paste0(
url_travel,
google_api_param,
google_api
)
url <- utils::URLencode(url, repeated = FALSE, reserved = FALSE)
message(paste0("Sent this request: ", url))
obj <- jsonlite::fromJSON(url)
if(obj$status != "OK" & any(grepl("error", names(obj))))
stop(obj[grepl("error", names(obj))], call. = FALSE)
if(grepl(pattern = "ZERO_RESULTS", obj$rows$elements[[1]])[1])
if (obj$status != "OK" & any(grepl("error", names(obj)))) {
stop(obj[grepl("error", names(obj))], call. = FALSE)
}
if (grepl(pattern = "ZERO_RESULTS", obj$rows$elements[[1]])[1]) {
stop("No results for this request (e.g. due to lack of support for this mode between the from and to locations)", call. = FALSE)
}

# some of cols are data.frames, e.g.
# lapply(obj$rows$elements[[1]], class)
# obj$rows$elements[[1]][1]
# obj$rows$elements[[1]][1]$distance$value
distances = lapply(obj$rows$elements,
function(x) x[1]$distance$value)
distances = unlist(distances)
duration = lapply(obj$rows$elements,
function(x) x[2]$duration$value)
duration = unlist(duration)
currency = NA
fare = NA
if(mode == "transit" & !is.null(obj$rows$elements[[1]]$fare[1])){
currency = lapply(obj$rows$elements,
function(x) x$fare$currency)
currency = unlist(currency)
fare = lapply(obj$rows$elements,
function(x) x$fare$value)
fare = unlist(fare)
distances <- lapply(
obj$rows$elements,
function(x) x[1]$distance$value
)
distances <- unlist(distances)
duration <- lapply(
obj$rows$elements,
function(x) x[2]$duration$value
)
duration <- unlist(duration)
currency <- NA
fare <- NA
if (mode == "transit" & !is.null(obj$rows$elements[[1]]$fare[1])) {
currency <- lapply(
obj$rows$elements,
function(x) x$fare$currency
)
currency <- unlist(currency)
fare <- lapply(
obj$rows$elements,
function(x) x$fare$value
)
fare <- unlist(fare)
}
# is_ok = lapply(obj$rows$elements,
# function(x) x$status)
from_addresses = rep(obj$origin_addresses, each = length(obj$origin_addresses))
to_addresses = rep(obj$destination_addresses, length(obj$origin_addresses))
res_df = data.frame(from_addresses, to_addresses, distances, duration, currency, fare)
from_addresses <- rep(obj$origin_addresses, each = length(obj$origin_addresses))
to_addresses <- rep(obj$destination_addresses, length(obj$origin_addresses))
res_df <- data.frame(from_addresses, to_addresses, distances, duration, currency, fare)
res_df$from_addresses <- as.character(res_df$from_addresses)
res_df$to_addresses <- as.character(res_df$to_addresses)
return(res_df)
}





116 R/gtfs.R
View file
@@ -14,16 +14,17 @@
#' # download.file(u, f)
#' gtfs <- gtfs2sldf(gtfszip = f)
#' plot(gtfs, col = gtfs$route_long_name)
#' plot(gtfs[gtfs$route_long_name == "Central Campus",])
#' plot(gtfs[gtfs$route_long_name == "Central Campus", ])
#' \dontrun{
#' # An example of a larger gtfs feed
#' download.file("http://www.yrt.ca/google/google_transit.zip",
#' paste0(tempdir(),"/gtfsfeed.zip"))
#' yrtgtfs <- gtfs2sldf(paste0(tempdir(),"/gtfsfeed.zip"))
#' sp::plot(yrtgtfs,col=paste0("#",yrtgtfs$route_color))
#' download.file(
#' "http://www.yrt.ca/google/google_transit.zip",
#' paste0(tempdir(), "/gtfsfeed.zip")
#' )
#' yrtgtfs <- gtfs2sldf(paste0(tempdir(), "/gtfsfeed.zip"))
#' sp::plot(yrtgtfs, col = paste0("#", yrtgtfs$route_color))
#' }
gtfs2sldf <- function(gtfszip = "") {

if (gtfszip == "") {
stop("Zip file required")
}
@@ -33,17 +34,17 @@ gtfs2sldf <- function(gtfszip = "") {

gtfsfiles <- unzip(gtfszip, exdir = tempdir())

gtfstrips <- read.csv(paste0(tempdir(),"/trips.txt"))
if (all(charToRaw(substr(colnames(gtfstrips)[1],1,3)) == c(as.raw(239),as.raw(46),as.raw(46)))) {
gtfstrips <- read.csv(paste0(tempdir(),"/trips.txt"),fileEncoding="UTF-8-BOM")
gtfsroutes <- read.csv(paste0(tempdir(),"/routes.txt"),fileEncoding="UTF-8-BOM")
gtfsagency <- read.csv(paste0(tempdir(),"/agency.txt"),fileEncoding="UTF-8-BOM")
gtfsshape <- read.csv(paste0(tempdir(),"/shapes.txt"),fileEncoding="UTF-8-BOM")
gtfstrips <- read.csv(paste0(tempdir(), "/trips.txt"))
if (all(charToRaw(substr(colnames(gtfstrips)[1], 1, 3)) == c(as.raw(239), as.raw(46), as.raw(46)))) {
gtfstrips <- read.csv(paste0(tempdir(), "/trips.txt"), fileEncoding = "UTF-8-BOM")
gtfsroutes <- read.csv(paste0(tempdir(), "/routes.txt"), fileEncoding = "UTF-8-BOM")
gtfsagency <- read.csv(paste0(tempdir(), "/agency.txt"), fileEncoding = "UTF-8-BOM")
gtfsshape <- read.csv(paste0(tempdir(), "/shapes.txt"), fileEncoding = "UTF-8-BOM")
}
else {
gtfsroutes <- read.csv(paste0(tempdir(),"/routes.txt"))
gtfsagency <- read.csv(paste0(tempdir(),"/agency.txt"))
gtfsshape <- read.csv(paste0(tempdir(),"/shapes.txt"))
gtfsroutes <- read.csv(paste0(tempdir(), "/routes.txt"))
gtfsagency <- read.csv(paste0(tempdir(), "/agency.txt"))
gtfsshape <- read.csv(paste0(tempdir(), "/shapes.txt"))
}

if (nrow(gtfsshape) == 0) {
@@ -53,46 +54,49 @@ gtfs2sldf <- function(gtfszip = "") {
unlink(gtfsfiles)

gtfslines <- sp::SpatialLinesDataFrame((gtfsshape %>%
dplyr::group_by_(~shape_id) %>%
dplyr::arrange_(~shape_pt_sequence) %>%
dplyr::do_(
gtfsline = "sp::Lines(sp::Line(as.matrix(.[,c('shape_pt_lon','shape_pt_lat')])),unique(.$shape_id))"
) %>%
dplyr::ungroup() %>%
dplyr::do_(
gtfsline = "sp::SpatialLines(.[[2]],
proj4string = sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))"))[[1]][[1]],
data=gtfstrips %>%
dplyr::inner_join(
gtfsroutes
) %>%
dplyr::distinct_(
~route_id,
~shape_id,
~route_short_name,
~route_long_name,
~route_desc,
~route_type,
~route_color,
~route_text_color,
~agency_id
) %>%
dplyr::select_(~route_id,
~shape_id,
~route_short_name,
~route_long_name,
~route_desc,
~route_type,
~route_color,
~route_text_color,
~agency_id) %>%
dplyr::inner_join(
gtfsagency
) %>%
dplyr::do_(
"`rownames<-`(.,.$shape_id)"
))
rm(gtfstrips,gtfsshape,gtfsagency,gtfsroutes)
dplyr::group_by_(~shape_id) %>%
dplyr::arrange_(~shape_pt_sequence) %>%
dplyr::do_(
gtfsline = "sp::Lines(sp::Line(as.matrix(.[,c('shape_pt_lon','shape_pt_lat')])),unique(.$shape_id))"
) %>%
dplyr::ungroup() %>%
dplyr::do_(
gtfsline = "sp::SpatialLines(.[[2]],
proj4string = sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))"
))[[1]][[1]],
data = gtfstrips %>%
dplyr::inner_join(
gtfsroutes
) %>%
dplyr::distinct_(
~route_id,
~shape_id,
~route_short_name,
~route_long_name,
~route_desc,
~route_type,
~route_color,
~route_text_color,
~agency_id
) %>%
dplyr::select_(
~route_id,
~shape_id,
~route_short_name,
~route_long_name,
~route_desc,
~route_type,
~route_color,
~route_text_color,
~agency_id
) %>%
dplyr::inner_join(
gtfsagency
) %>%
dplyr::do_(
"`rownames<-`(.,.$shape_id)"
)
)
rm(gtfstrips, gtfsshape, gtfsagency, gtfsroutes)
return(gtfslines)
}

24 R/line_match.R
View file
@@ -8,30 +8,30 @@
#' @param return_sp Should the function return a spatial result (FALSE by default)
#' @export
#' @examples
#' x1 = 2:4
#' x2 = 3:5
#' x1 <- 2:4
#' x2 <- 3:5
#' match(x1, x2) # how the base function works
#' l1 = flowlines[2:4,]
#' l2 = routes_fast[3:5,]
#' (lmatches = line_match(l1, l2)) # how the stplanr version works
#' l2matched = l2[lmatches[!is.na(lmatches)],]
#' l1 <- flowlines[2:4, ]
#' l2 <- routes_fast[3:5, ]
#' (lmatches <- line_match(l1, l2)) # how the stplanr version works
#' l2matched <- l2[lmatches[!is.na(lmatches)], ]
#' plot(l1)
#' plot(l2, add = TRUE)
#' plot(l2matched, add = TRUE, col = "red") # showing matched routes
#' l2matched2 = line_match(l1, l2, return_sp = TRUE)
#' l2matched2 <- line_match(l1, l2, return_sp = TRUE)
#' identical(l2matched, l2matched2)
#' # decreasing the match likelihood via the threshold
#' line_match(l1, l2, threshold = 0.003)
line_match <- function(l1, l2, threshold = 0.01, return_sp = FALSE){
line_match <- function(l1, l2, threshold = 0.01, return_sp = FALSE) {
dist_mat <- rgeos::gDistance(l1, l2, byid = TRUE, hausdorff = TRUE)
closest <- apply(dist_mat, 2, which.min)
closest_values <- apply(dist_mat, 2, min)
closest[closest_values > threshold] <- NA
sel_na <- is.na(closest)
if(return_sp) {
l2matched <- l2[closest[!sel_na],]
if (return_sp) {
l2matched <- l2[closest[!sel_na], ]
match_num <- names(closest)[which(!sel_na)]
if(is(l2, "SpatialLinesDataFrame")) {
if (is(l2, "SpatialLinesDataFrame")) {
l2matched$match <- match_num
} else {
l2matched <- sp::SpatialLinesDataFrame(l2matched, data = data.frame(match_num), match.ID = FALSE)
@@ -40,4 +40,4 @@ line_match <- function(l1, l2, threshold = 0.01, return_sp = FALSE){
} else {
return(unname(closest))
}
}
}
91 R/line_sample.R
View file
@@ -5,34 +5,35 @@
#' @param weights Relative probabilities of samples on lines
#' @export
#' @examples
#' n = 10
#' l_lengths = 1:5
#' weights = 9:5
#' (res = n_sample_length(n, l_lengths, weights))
#' n <- 10
#' l_lengths <- 1:5
#' weights <- 9:5
#' (res <- n_sample_length(n, l_lengths, weights))
#' sum(res)
#' n = 100
#' l_lengths = c(12, 22, 15, 14)
#' weights = c(38, 10, 44, 34)
#' (res = n_sample_length(n, l_lengths, weights))
#' n <- 100
#' l_lengths <- c(12, 22, 15, 14)
#' weights <- c(38, 10, 44, 34)
#' (res <- n_sample_length(n, l_lengths, weights))
#' sum(res)
#' # more examples:
#' n_sample_length(5, 1:5, c(0.1, 0.9, 0, 0, 0))
#' n_sample_length(5, 1:5, c(0.5, 0.3, 0.1, 0, 0))
#' l = flowlines[2:6,]
#' l_lengths = line_length(l)
#' n = n_sample_length(10, l_lengths, weights = l$All)
#' l <- flowlines[2:6, ]
#' l_lengths <- line_length(l)
#' n <- n_sample_length(10, l_lengths, weights = l$All)
n_sample_length <- function(n, l_lengths, weights) {
# generate length-adjusted weights equal to 1
l_lengths_rel = l_lengths * weights / (sum(l_lengths * weights))
n_vec = round(n * l_lengths_rel)
if(sum(n_vec) != n) {
n_diff = n - sum(n_vec)
probs = l_lengths_rel - (n_vec / n) # how much less of each is needed
if(n_diff < 0)
probs = probs * -1
probs[probs < 0] = 0
sel = sample(length(l_lengths), size = abs(n_diff), prob = probs)
n_vec[sel] = n_vec[sel] + sign(n_diff)
l_lengths_rel <- l_lengths * weights / (sum(l_lengths * weights))
n_vec <- round(n * l_lengths_rel)
if (sum(n_vec) != n) {
n_diff <- n - sum(n_vec)
probs <- l_lengths_rel - (n_vec / n) # how much less of each is needed
if (n_diff < 0) {
probs <- probs * -1
}
probs[probs < 0] <- 0
sel <- sample(length(l_lengths), size = abs(n_diff), prob = probs)
n_vec[sel] <- n_vec[sel] + sign(n_diff)
}
n_vec
}
@@ -44,35 +45,37 @@ n_sample_length <- function(n, l_lengths, weights) {
#' @param weights The relative probabilities of lines being samples
#' @export
#' @examples
#' l = flowlines[2:5,]
#' n = 100
#' l_lengths = line_length(l)
#' weights = l$All
#' p = line_sample(l, 50, weights)
#' l <- flowlines[2:5, ]
#' n <- 100
#' l_lengths <- line_length(l)
#' weights <- l$All
#' p <- line_sample(l, 50, weights)
#' plot(p)
#' p = line_sample(l, 50, weights = 1:length(l))
#' p <- line_sample(l, 50, weights = 1:length(l))
#' plot(p)
line_sample <- function(l, n, weights) {
not_projected = !sp::is.projected(l)
if(not_projected) {
crs_orig = sp::proj4string(l)
if(is.na(crs_orig))
crs_orig = sp::CRS("+init=epsg:4326")
crs_new = crs_select_aeq(l)
l = sp::spTransform(l, CRSobj = crs_new)
not_projected <- !sp::is.projected(l)
if (not_projected) {
crs_orig <- sp::proj4string(l)
if (is.na(crs_orig)) {
crs_orig <- sp::CRS("+init=epsg:4326")
}
crs_new <- crs_select_aeq(l)
l <- sp::spTransform(l, CRSobj = crs_new)
}
lsf = sf::st_as_sf(l, "SpatialLinesDataFrame")
l_lengths = sf::st_length(lsf)
l_lengths = as.vector(l_lengths) # convert to numeric vector
n_vec = n_sample_length(n, l_lengths, weights = weights)
psf = sf::st_line_sample(lsf, n = n_vec)
lsf <- sf::st_as_sf(l, "SpatialLinesDataFrame")
l_lengths <- sf::st_length(lsf)
l_lengths <- as.vector(l_lengths) # convert to numeric vector
n_vec <- n_sample_length(n, l_lengths, weights = weights)
psf <- sf::st_line_sample(lsf, n = n_vec)

# aim: group point collection into points - to update if possible
# class(psf) # its a MULIPOINT
psf_point = as(psf, "Spatial")
psp = as(psf, "Spatial")
psp = sp::SpatialPoints(matrix(sp::coordinates(psp), ncol = 2), proj4string = sp::CRS(proj4string(psp)))
if(not_projected)
psp = sp::spTransform(psp, crs_orig)
psf_point <- as(psf, "Spatial")
psp <- as(psf, "Spatial")
psp <- sp::SpatialPoints(matrix(sp::coordinates(psp), ncol = 2), proj4string = sp::CRS(proj4string(psp)))
if (not_projected) {
psp <- sp::spTransform(psp, crs_orig)
}
return(psp)
}
110 R/linefuns.R
View file
@@ -13,15 +13,15 @@
#' @examples
#' n_vertices(routes_fast)
#' n_vertices(routes_fast_sf)
n_vertices <- function(l){
n_vertices <- function(l) {
UseMethod("n_vertices")
}
#' @export
n_vertices.Spatial <- function(l){
n_vertices.Spatial <- function(l) {
sapply(l@lines, function(x) nrow(x@Lines[[1]]@coords))
}
#' @export
n_vertices.sf <- function(l){
n_vertices.sf <- function(l) {
geoms <- sf::st_coordinates(l)
L1 <- rlang::quo(L1)
geoms %>%
@@ -49,8 +49,8 @@ n_vertices.sf <- function(l){
#' nrow(flowlines)
#' sum(islp)
#' # Remove invisible 'linepoints'
#' nrow(flowlines[!islp,])
is_linepoint <- function(l){
#' nrow(flowlines[!islp, ])
is_linepoint <- function(l) {
nverts <- n_vertices(l)
sel <- nverts <= 2
ldf <- line2df(l)
@@ -80,12 +80,11 @@ line_bearing <- function(l, bidirectional = FALSE) {
}
#' @export
line_bearing.Spatial <- function(l, bidirectional = FALSE) {

ldf <- line2df(l)
bearing <- geosphere::bearing(as.matrix(ldf[, c("fx", "fy")]), as.matrix(ldf[,c("tx", "ty")]))
if(bidirectional) {
bearing[bearing > 90] <- bearing[bearing > 90] - 180
bearing[bearing < -90] <- bearing[bearing < -90] + 180
bearing <- geosphere::bearing(as.matrix(ldf[, c("fx", "fy")]), as.matrix(ldf[, c("tx", "ty")]))
if (bidirectional) {
bearing[bearing > 90] <- bearing[bearing > 90] - 180
bearing[bearing < -90] <- bearing[bearing < -90] + 180
}
bearing
}
@@ -113,35 +112,36 @@ line_bearing.sf <- function(l, bidirectional = FALSE) {
#' @examples
#' data(flowlines)
#' # Find all routes going North-South
#' a = angle_diff(flowlines, angle = 0, bidirectional = TRUE, absolute = TRUE)
#' a <- angle_diff(flowlines, angle = 0, bidirectional = TRUE, absolute = TRUE)
#' plot(flowlines)
#' plot(flowlines[a < 15,], add = TRUE, lwd = 3, col = "red")
#' plot(flowlines[a < 15, ], add = TRUE, lwd = 3, col = "red")
#' # East-West
#' plot(flowlines[a > 75,], add = TRUE, lwd = 3, col = "green")
#' plot(flowlines[a > 75, ], add = TRUE, lwd = 3, col = "green")
#' angle_diff(flowlines_sf[2, ], angle = 0)
angle_diff <- function(l, angle, bidirectional = FALSE, absolute = TRUE) {
UseMethod("angle_diff")
}
#' @export
angle_diff.Spatial <- function(l, angle, bidirectional = FALSE, absolute = TRUE){
if(is(object = l, "Spatial")){
line_angles = line_bearing(l)
angle_diff.Spatial <- function(l, angle, bidirectional = FALSE, absolute = TRUE) {
if (is(object = l, "Spatial")) {
line_angles <- line_bearing(l)
} else {
line_angles = l
line_angles <- l
}
angle_diff <- angle - line_angles
angle_diff[angle_diff <= -180] <- angle_diff[angle_diff <= -180] + 180
angle_diff[angle_diff >= 180] <- angle_diff[angle_diff >= 180] - 180
if (bidirectional) {
angle_diff[angle_diff <= -90] <- 180 + angle_diff[angle_diff <= -90]
angle_diff[angle_diff >= 90] <- 180 - angle_diff[angle_diff >= 90]
}
angle_diff = angle - line_angles
angle_diff[angle_diff <= -180] = angle_diff[angle_diff <= -180] + 180
angle_diff[angle_diff >= 180] = angle_diff[angle_diff >= 180] - 180
if(bidirectional){
angle_diff[angle_diff <= -90] = 180 + angle_diff[angle_diff <= -90]
angle_diff[angle_diff >= 90] = 180 - angle_diff[angle_diff >= 90]
if (absolute) {
angle_diff <- abs(angle_diff)
}
if(absolute)
angle_diff = abs(angle_diff)
angle_diff
}
#' @export
angle_diff.sf <- function(l, angle, bidirectional = FALSE, absolute = TRUE){
angle_diff.sf <- function(l, angle, bidirectional = FALSE, absolute = TRUE) {
l_sp <- as(l, "Spatial")
angle_diff.Spatial(l_sp, angle, bidirectional = FALSE, absolute = TRUE)
}
@@ -153,7 +153,7 @@ angle_diff.sf <- function(l, angle, bidirectional = FALSE, absolute = TRUE){
#' @export
#' @examples
#' data(routes_fast)
#' line_midpoint(routes_fast[2:5,])
#' line_midpoint(routes_fast[2:5, ])
line_midpoint <- function(l) {
UseMethod("line_midpoint")
}
@@ -171,7 +171,7 @@ line_midpoint.sf <- function(l) {
#' @inheritParams line2df
#' @param byid Logical determining whether the length is returned per object (default is true)
#' @export
line_length <- function(l, byid = TRUE){
line_length <- function(l, byid = TRUE) {
gprojected(l, rgeos::gLength, byid = byid)
}

@@ -182,40 +182,40 @@ line_length <- function(l, byid = TRUE){
#' @export
#' @examples
#' data(routes_fast)
#' l = routes_fast[2, ]
#' l <- routes_fast[2, ]
#' library(sp)
#' l_seg2 = line_segment(l = l, n_segments = 2)
#' l_seg2 <- line_segment(l = l, n_segments = 2)
#' plot(l_seg2, col = l_seg2$group, lwd = 50)
line_segment <- function(l, n_segments, segment_length = NA){
if(!is.na(segment_length)){
l_length = line_length(l)
n_segments = round(l_length / segment_length)
line_segment <- function(l, n_segments, segment_length = NA) {
if (!is.na(segment_length)) {
l_length <- line_length(l)
n_segments <- round(l_length / segment_length)
}
if(n_segments == 2){
pseg = line_midpoint(l)
if (n_segments == 2) {
pseg <- line_midpoint(l)
} else {
pseg = sp::spsample(x = l, n = n_segments - 1, type = "regular")
pseg <- sp::spsample(x = l, n = n_segments - 1, type = "regular")
}
l_geom = raster::geom(l)
l_coords = l_geom[, c("x", "y")]
knn_res = nabor::knn(data = l_coords, query = sp::coordinates(pseg), k = 1)
sel_nearest = c(knn_res$nn.idx)
for(i in 1:(length(sel_nearest) + 1)){
ids = c(1, sel_nearest, nrow(l))
if(i == 1){
l_seg = points2line(l_coords[ids[i]:ids[(i + 1)],])
sp::spChFIDs(l) = i
} else if(i == length(sel_nearest) + 1){
l_temp = points2line(l_coords[ids[i]:nrow(l_coords),])
sp::spChFIDs(l_temp) = i
l_seg = raster::bind(l_seg, l_temp)
l_geom <- raster::geom(l)
l_coords <- l_geom[, c("x", "y")]
knn_res <- nabor::knn(data = l_coords, query = sp::coordinates(pseg), k = 1)
sel_nearest <- c(knn_res$nn.idx)
for (i in 1:(length(sel_nearest) + 1)) {
ids <- c(1, sel_nearest, nrow(l))
if (i == 1) {
l_seg <- points2line(l_coords[ids[i]:ids[(i + 1)], ])
sp::spChFIDs(l) <- i
} else if (i == length(sel_nearest) + 1) {
l_temp <- points2line(l_coords[ids[i]:nrow(l_coords), ])
sp::spChFIDs(l_temp) <- i
l_seg <- raster::bind(l_seg, l_temp)
} else {
l_temp = points2line(l_coords[ids[i]:ids[(i + 1)],])
sp::spChFIDs(l_temp) = i
l_seg = raster::bind(l_seg, l_temp)
l_temp <- points2line(l_coords[ids[i]:ids[(i + 1)], ])
sp::spChFIDs(l_temp) <- i
l_seg <- raster::bind(l_seg, l_temp)
}
}
l_seg = sp::SpatialLinesDataFrame(l_seg, data.frame(group = 1:i))
raster::crs(l_seg) = raster::crs(l)
l_seg <- sp::SpatialLinesDataFrame(l_seg, data.frame(group = 1:i))
raster::crs(l_seg) <- raster::crs(l)
l_seg
}
239 R/load-stats19.R
View file

Large diffs are not rendered by default.

69 R/loadABS.R
View file
@@ -23,22 +23,23 @@
#' @export
#' @examples
#' data_dir <- system.file("extdata", package = "stplanr")
#' t1 <- read_table_builder(file.path(data_dir, 'SA1Population.csv'))
#' t2 <- read_table_builder(file.path(data_dir, 'SA1Population.xlsx'),
#' filetype = 'xlsx', sheet = 1, removeTotal = TRUE)
#' sa1pop <- read.csv(file.path(data_dir, 'SA1Population.csv'), header=FALSE)
#' t1 <- read_table_builder(file.path(data_dir, "SA1Population.csv"))
#' t2 <- read_table_builder(file.path(data_dir, "SA1Population.xlsx"),
#' filetype = "xlsx", sheet = 1, removeTotal = TRUE
#' )
#' sa1pop <- read.csv(file.path(data_dir, "SA1Population.csv"), header = FALSE)
#' t3 <- read_table_builder(sa1pop)
read_table_builder <- function(dataset, filetype="csv",sheet=1,removeTotal=TRUE) {
read_table_builder <- function(dataset, filetype = "csv", sheet = 1, removeTotal = TRUE) {
if (missing(dataset)) {
stop("Dataset is missing")
}
if (is.data.frame(dataset)) {
tbfile <- dataset
} else if (is.character(dataset)) {
if (filetype=="xlsx") {
tbfile <- openxlsx::readWorkbook(dataset,sheet = sheet, colNames = FALSE)
if (filetype == "xlsx") {
tbfile <- openxlsx::readWorkbook(dataset, sheet = sheet, colNames = FALSE)
} else {
tbfile <- read.csv(dataset,header=FALSE)
tbfile <- read.csv(dataset, header = FALSE)
}
} else {
stop("Dataset not data.frame or character string")
@@ -47,52 +48,54 @@ read_table_builder <- function(dataset, filetype="csv",sheet=1,removeTotal=TRUE)
stop("File could not be loaded")
} else {
if (filetype == "xlsx" | filetype == "legacycsv") {
tbfile[tbfile == ''] <- NA
tbfile <- tbfile[,which(!(colSums(is.na(tbfile)) == nrow(tbfile)))]
if (is.na(tbfile[which(rowSums(is.na(tbfile[,2:ncol(tbfile)])) == min(rowSums(is.na(tbfile[,2:ncol(tbfile)])))),][1,1]) == TRUE) {
tbfile[,1] <- NULL
tbfile[tbfile == ""] <- NA
tbfile <- tbfile[, which(!(colSums(is.na(tbfile)) == nrow(tbfile)))]
if (is.na(tbfile[which(rowSums(is.na(tbfile[, 2:ncol(tbfile)])) == min(rowSums(is.na(tbfile[, 2:ncol(tbfile)])))), ][1, 1]) == TRUE) {
tbfile[, 1] <- NULL
}
else {
tbfile <- tbfile[which(rowSums(is.na(tbfile)) < (ncol(tbfile)-1)),]
tbfile <- tbfile[which(rowSums(is.na(tbfile)) < (ncol(tbfile) - 1)), ]
}
tbfile <- tbfile[which(rowSums(is.na(tbfile)) != ncol(tbfile)),]
valuecols <- which(!is.na(tbfile[1,]))
tbfile <- tbfile[which(rowSums(is.na(tbfile)) != ncol(tbfile)), ]
valuecols <- which(!is.na(tbfile[1, ]))
valuecols <- valuecols[which(valuecols > 1)]
valuecols <- valuecols[!valuecols %in% which(!is.na(tbfile[2,]))]
colnames(tbfile) <- c(as.character(unlist(unname(tbfile[2,which(!is.na(tbfile[2,]))]))),as.character(unlist(unname(tbfile[1,valuecols]))))
tbfile <- tbfile[3:nrow(tbfile),]
tbfile <- tbfile[which(rowSums(is.na(tbfile)) != ncol(tbfile)-1),]
valuecols <- valuecols[!valuecols %in% which(!is.na(tbfile[2, ]))]
colnames(tbfile) <- c(as.character(unlist(unname(tbfile[2, which(!is.na(tbfile[2, ]))]))), as.character(unlist(unname(tbfile[1, valuecols]))))
tbfile <- tbfile[3:nrow(tbfile), ]
tbfile <- tbfile[which(rowSums(is.na(tbfile)) != ncol(tbfile) - 1), ]
if (length(valuecols) > 1) {
tbfile <- tbfile[which(!rowSums(is.na(tbfile[,valuecols])) == length(valuecols)),]
tbfile <- tbfile[which(!rowSums(is.na(tbfile[, valuecols])) == length(valuecols)), ]
}
else {
tbfile <- tbfile[which(is.na(tbfile[,valuecols]) != TRUE),]
tbfile <- tbfile[which(is.na(tbfile[, valuecols]) != TRUE), ]
}
i <- 1
while (sum(is.na(tbfile[,i])) != 0) {
tbfile[,i] <- rep(
unique(tbfile[which(is.na(tbfile[,i])==FALSE),i]),
each=nrow(tbfile)/length(tbfile[which(is.na(tbfile[,i])==FALSE),i]),
times=length(tbfile[which(is.na(tbfile[,i])==FALSE),i])/length(unique(tbfile[which(is.na(tbfile[,i])==FALSE),i]))
while (sum(is.na(tbfile[, i])) != 0) {
tbfile[, i] <- rep(
unique(tbfile[which(is.na(tbfile[, i]) == FALSE), i]),
each = nrow(tbfile) / length(tbfile[which(is.na(tbfile[, i]) == FALSE), i]),
times = length(tbfile[which(is.na(tbfile[, i]) == FALSE), i]) / length(unique(tbfile[which(is.na(tbfile[, i]) == FALSE), i]))
)
i <- i + 1
}
if (removeTotal == TRUE) {
tbfile <- tbfile[,which(colnames(tbfile) != "Total")]
tbfile <- tbfile[which(tbfile[,1] != "Total"),]
tbfile <- tbfile[, which(colnames(tbfile) != "Total")]
tbfile <- tbfile[which(tbfile[, 1] != "Total"), ]
}
tbfile[valuecols[which(valuecols <= ncol(tbfile))]] <- sapply(tbfile[valuecols[which(valuecols <= ncol(tbfile))]],function(x){as.numeric(as.character(x))})
tbfile[valuecols[which(valuecols <= ncol(tbfile))]] <- sapply(tbfile[valuecols[which(valuecols <= ncol(tbfile))]], function(x) {
as.numeric(as.character(x))
})
row.names(tbfile) <- NULL
} else {
colnamevals <- c(as.character(unname(unlist(tbfile[(min(which(is.na(tbfile[,ncol(tbfile)])==FALSE))-1),1:(ncol(tbfile)-1)]))),'value')
tbfile <- tbfile[which(is.na(tbfile[,ncol(tbfile)])==FALSE),]
colnamevals <- c(as.character(unname(unlist(tbfile[(min(which(is.na(tbfile[, ncol(tbfile)]) == FALSE)) - 1), 1:(ncol(tbfile) - 1)]))), "value")
tbfile <- tbfile[which(is.na(tbfile[, ncol(tbfile)]) == FALSE), ]
colnames(tbfile) <- colnamevals
if (removeTotal == TRUE) {
tbfile <- tbfile[apply(tbfile, 1, function(x) all(x != 'Total')),]
tbfile <- tbfile[apply(tbfile, 1, function(x) all(x != "Total")), ]
}
row.names(tbfile) <- NULL
tbfile$value <- as.numeric(as.character(tbfile$value))
}
}
return(tbfile)
}
}
289 R/od-funs.R
View file

Large diffs are not rendered by default.

24 R/od_coords.R
View file
@@ -17,42 +17,40 @@
#' od_coords(flowlines[1:3, ])
#' od_coords(flowlines_sf[1:3, ])
od_coords <- function(from = NULL, to = NULL, l = NULL) {

if(is(object = from, class2 = "sf")) {
if (is(object = from, class2 = "sf")) {
is_sf_line <- all(sf::st_geometry_type(from) == "LINESTRING")
} else {
is_sf_line <- FALSE
}

if(is_sf_line | any(grepl(pattern = "Line", x = class(from)))) {
if (is_sf_line | any(grepl(pattern = "Line", x = class(from)))) {
l <- from
}

if(!is.null(l)) {
if (!is.null(l)) {
coord_matrix <- line2df(l) %>%
dplyr::select("fx", "fy", "tx", "ty")
}

else {
# Convert sp object to lat/lon vector
if(is(object = from, "Spatial")) from <- sp::coordinates(from)
if(is(object = to, "Spatial")) to <- sp::coordinates(to)
if (is(object = from, "Spatial")) from <- sp::coordinates(from)
if (is(object = to, "Spatial")) to <- sp::coordinates(to)

# sf objects
if(is(object = from, "sf")) from <- sf::st_coordinates(from)
if(is(object = to, "sf")) to <- sf::st_coordinates(to)
if (is(object = from, "sf")) from <- sf::st_coordinates(from)
if (is(object = to, "sf")) to <- sf::st_coordinates(to)

# Convert character strings to lon/lat if needs be
if(is.character(from)) from <- matrix(geo_code(from), ncol = 2)
if(is.character(to)) to <- matrix(geo_code(to), ncol = 2)
if(is.vector(from) & is.vector(to)) {
coord_matrix = matrix(c(from, to), ncol = 4)
if (is.character(from)) from <- matrix(geo_code(from), ncol = 2)
if (is.character(to)) to <- matrix(geo_code(to), ncol = 2)
if (is.vector(from) & is.vector(to)) {
coord_matrix <- matrix(c(from, to), ncol = 4)
} else {
coord_matrix <- cbind(from, to)
}
colnames(coord_matrix) <- c("fx", "fy", "tx", "ty")
}

coord_matrix

}
60 R/oneway.R
View file
@@ -33,28 +33,27 @@ onewayid <- function(x, attrib, id1 = names(x)[1], id2 = names(x)[2],
#' for users and are difficult to plot.
#' @examples
#' data(flow)
#' flow_oneway = onewayid(flow, attrib = 3)
#' flow_oneway <- onewayid(flow, attrib = 3)
#' nrow(flow_oneway) < nrow(flow) # result has fewer rows
#' sum(flow$All) == sum(flow_oneway$All) # but the same total flow
#' # using names instead of index for attribute
#' onewayid(flow, attrib = "All")
#' # using many attributes to aggregate
#' attrib = which(vapply(flow, is.numeric, TRUE))
#' flow_oneway = onewayid(flow, attrib = attrib)
#' attrib <- which(vapply(flow, is.numeric, TRUE))
#' flow_oneway <- onewayid(flow, attrib = attrib)
#' colSums(flow_oneway[attrib]) == colSums(flow[attrib]) # test if the colSums are equal
#' # Demonstrate the results from onewayid and onewaygeo are identical
#' flow_oneway_geo = onewaygeo(flowlines, attrib = attrib)
#' flow_oneway_geo <- onewaygeo(flowlines, attrib = attrib)
#' plot(flow_oneway$All, flow_oneway_geo$All)
#' onewayid(flowlines_sf, "all")
#' @export
onewayid.data.frame <- function(x, attrib, id1 = names(x)[1], id2 = names(x)[2],
stplanr.key = od_id_order(x, id1, id2)){

if(is.numeric(attrib)){
attrib_names = names(x)[attrib]
stplanr.key = od_id_order(x, id1, id2)) {
if (is.numeric(attrib)) {
attrib_names <- names(x)[attrib]
} else {
attrib_names = attrib
attrib = which(names(x) %in% attrib)
attrib_names <- attrib
attrib <- which(names(x) %in% attrib)
}

x <- dplyr::bind_cols(x, stplanr.key)
@@ -67,27 +66,26 @@ onewayid.data.frame <- function(x, attrib, id1 = names(x)[1], id2 = names(x)[2],
dplyr::summarise(is_two_way = dplyr::last(.data$is_two_way)) %>%
dplyr::select(-stplanr.key)

x_oneway_character = x %>%
x_oneway_character <- x %>%
dplyr::transmute(
id1 = stringr::str_split(.data$stplanr.key, " ", simplify = TRUE)[, 1],
id2 = stringr::str_split(.data$stplanr.key, " ", simplify = TRUE)[, 2],
stplanr.key = .data$stplanr.key
) %>%
dplyr::group_by(stplanr.key) %>%
dplyr::summarise(id1 =dplyr::first(id1), id2 =dplyr::first(id2)) %>%
dplyr::summarise(id1 = dplyr::first(id1), id2 = dplyr::first(id2)) %>%
dplyr::select(-stplanr.key)

x_oneway = dplyr::bind_cols(
x_oneway <- dplyr::bind_cols(
x_oneway_character,
x_oneway_numeric,
x_oneway_binary
)

x_oneway$stplanr.key <- NULL
names(x_oneway)[1:2] = c(id1, id2)
names(x_oneway)[1:2] <- c(id1, id2)

return(x_oneway)

}

#' @name onewayid
@@ -100,43 +98,43 @@ onewayid.data.frame <- function(x, attrib, id1 = names(x)[1], id2 = names(x)[2],
#' sum(fo$All) == sum(flowlines$All)
#' # test results for one line
#' n <- 3
#' plot(fo[n,], lwd = 20, add = TRUE)
#' f_over_n <- rgeos::gEquals(fo[n,], flowlines["All"], byid = TRUE)
#' plot(fo[n, ], lwd = 20, add = TRUE)
#' f_over_n <- rgeos::gEquals(fo[n, ], flowlines["All"], byid = TRUE)
#' sum(flowlines$All[f_over_n]) == sum(fo$All[n]) # check aggregation worked
#' plot(flowlines[which(f_over_n)[1],], add = TRUE, col = "white", lwd = 10)
#' plot(flowlines[which(f_over_n)[2],], add = TRUE, lwd = 5)
#' plot(flowlines[which(f_over_n)[1], ], add = TRUE, col = "white", lwd = 10)
#' plot(flowlines[which(f_over_n)[2], ], add = TRUE, lwd = 5)
#' @export
onewayid.SpatialLines <- function(x, attrib, id1 = names(x)[1], id2 = names(x)[2],
stplanr.key = od_id_order(x, id1, id2)){

stplanr.key = od_id_order(x, id1, id2)) {
x_geom <- sp::SpatialLines(x@lines, proj4string = sp::CRS(proj4string(x)))
x <- x@data

x_oneway = onewayid(x, id1, id2, attrib = attrib)
x_oneway$stplanr.key = od_id_order(x_oneway[c(id1, id2)])$stplanr.key
x_oneway <- onewayid(x, id1, id2, attrib = attrib)
x_oneway$stplanr.key <- od_id_order(x_oneway[c(id1, id2)])$stplanr.key

if(length(x_geom) != nrow(x_oneway)) {
if (length(x_geom) != nrow(x_oneway)) {
id_old <- paste(x[[id1]], x[[id2]])
sel <- id_old %in% x_oneway$stplanr.key
x_geom <- x_geom[sel,]
x_geom <- x_geom[sel, ]
}

x_oneway <- sp::SpatialLinesDataFrame(sl = x_geom, data = x_oneway, match.ID = FALSE)

return(x_oneway)

}

#' Generate ordered ids of OD pairs so lowest is always first
#'
#' @inheritParams onewayid
#'
#' @examples
#' x = data.frame(id1 = c(1, 1, 2, 2, 3), id2 = c(1, 2, 3, 1, 4))
#' x <- data.frame(id1 = c(1, 1, 2, 2, 3), id2 = c(1, 2, 3, 1, 4))
#' od_id_order(x) # 4th line switches id1 and id2 so stplanr.key is in order
#' @export
od_id_order <- function(x, id1 = names(x)[1], id2 = names(x)[2]){
dplyr::transmute_(x, stplanr.id1 = as.name(id1),
stplanr.id2 = as.name(id2),
stplanr.key = ~paste(pmin(stplanr.id1, stplanr.id2), pmax(stplanr.id1, stplanr.id2)))
od_id_order <- function(x, id1 = names(x)[1], id2 = names(x)[2]) {
dplyr::transmute_(x,
stplanr.id1 = as.name(id1),
stplanr.id2 = as.name(id2),
stplanr.key = ~paste(pmin(stplanr.id1, stplanr.id2), pmax(stplanr.id1, stplanr.id2))
)
}
837 R/osrmfunctions.R
View file

Large diffs are not rendered by default.

202 R/overline.R
View file
@@ -7,22 +7,23 @@
#' @param g1 A spatial object
#' @param g2 A spatial object
#' @export
#' @examples \dontrun{
#' rnet <- overline(routes_fast[c(2, 3, 22),], attrib = "length")
#' @examples
#' \dontrun{
#' rnet <- overline(routes_fast[c(2, 3, 22), ], attrib = "length")
#' plot(rnet)
#' lines(routes_fast[22,], col = "red") # line without overlaps
#' islines(routes_fast[2,], routes_fast[3,])
#' islines(routes_fast[2,], routes_fast[22,])
#' lines(routes_fast[22, ], col = "red") # line without overlaps
#' islines(routes_fast[2, ], routes_fast[3, ])
#' islines(routes_fast[2, ], routes_fast[22, ])
#' # sf implementation
#' islines(routes_fast_sf[2,], routes_fast_sf[3,])
#' islines(routes_fast_sf[2,], routes_fast_sf[22,])
#' islines(routes_fast_sf[2, ], routes_fast_sf[3, ])
#' islines(routes_fast_sf[2, ], routes_fast_sf[22, ])
#' }
islines <- function(g1, g2) {
UseMethod("islines")
}
islines.Spatial <- function(g1, g2){
islines.Spatial <- function(g1, g2) {
## return TRUE if geometries intersect as lines, not points
inherits(rgeos::gIntersection(g1,g2), "SpatialLines")
inherits(rgeos::gIntersection(g1, g2), "SpatialLines")
}
islines.sf <- function(g1, g2) {
sf::st_geometry_type(sf::st_intersection(sf::st_geometry(g1), sf::st_geometry(g2))) == "MULTILINESTRING"
@@ -39,43 +40,41 @@ islines.sf <- function(g1, g2) {
#' @param buff_dist A number specifying the distance in meters of the buffer to be used to crop lines before running the operation. If the distance is zero (the default) touching but non-overlapping lines may be aggregated.
#' @export
#' @examples
#' sl <- routes_fast[2:4,]
#' sl <- routes_fast[2:4, ]
#' rsec <- gsection(sl)
#' rsec_buff <- gsection(sl, buff_dist = 1)
#' plot(sl[1], lwd = 9, col = 1:nrow(sl))
#' plot(rsec, col = 5 + (1:length(rsec)), add = TRUE, lwd = 3)
#' plot(rsec_buff, col = 5 + (1:length(rsec_buff)), add = TRUE, lwd = 3)
#' \dontrun{
#' sl <- routes_fast_sf[2:4,]
#' rsec <- gsection(sl)
#' rsec <- gsection(sl, buff_dist = 100) # 4 features: issue
#' sl <- routes_fast_sf[2:4, ]
#' rsec <- gsection(sl)
#' rsec <- gsection(sl, buff_dist = 100) # 4 features: issue
#' }
gsection <- function(sl, buff_dist = 0) {
UseMethod("gsection")
}
#' @export
gsection.Spatial <- function(sl, buff_dist = 0){
if(buff_dist > 0) {
sl = geo_toptail(sl, toptail_dist = buff_dist)
gsection.Spatial <- function(sl, buff_dist = 0) {
if (buff_dist > 0) {
sl <- geo_toptail(sl, toptail_dist = buff_dist)
}
overlapping = rgeos::gOverlaps(sl, byid = T)
overlapping <- rgeos::gOverlaps(sl, byid = T)
u <- rgeos::gUnion(sl, sl)
u_merged <- rgeos::gLineMerge(u)
sp::disaggregate(u_merged)
}
#' @export
gsection.sf <- function(sl, buff_dist = 0){

if(buff_dist > 0) {
sl = geo_toptail(sl, toptail_dist = buff_dist)
gsection.sf <- function(sl, buff_dist = 0) {
if (buff_dist > 0) {
sl <- geo_toptail(sl, toptail_dist = buff_dist)
}

u <- sf::st_union(sl)
u_merged <- sf::st_line_merge(u)
u_disag <- sf::st_cast(u_merged, "LINESTRING")

u_disag

}
#' Label SpatialLinesDataFrame objects
#'
@@ -89,10 +88,10 @@ gsection.sf <- function(sl, buff_dist = 0){
#' @author Barry Rowlingson
#'
#' @export
lineLabels <- function(sl, attrib){
lineLabels <- function(sl, attrib) {
text(sp::coordinates(
rgeos::gCentroid(sl, byid = TRUE)
), labels = sl[[attrib]])
), labels = sl[[attrib]])
}

#' Convert series of overlapping lines into a route network
@@ -117,28 +116,29 @@ lineLabels <- function(sl, attrib){
#' \url{http://gis.stackexchange.com}. See \url{http://gis.stackexchange.com/questions/139681/overlaying-lines-and-aggregating-their-values-for-overlapping-segments}.
#' @export
#' @examples
#' sl <- routes_fast[2:4,]
#' sl <- routes_fast[2:4, ]
#' rnet1 <- overline(sl = sl, attrib = "length")
#' rnet2 <- overline(sl = sl, attrib = "length", buff_dist = 1)
#' plot(rnet1, lwd = rnet1$length / mean(rnet1$length))
#' plot(rnet2, lwd = rnet2$length / mean(rnet2$length))
#' \dontrun{
#' routes_fast$group = rep(1:3, length.out = nrow(routes_fast))
#' rnet_grouped = overline(routes_fast, attrib = "length", byvars = "group", buff_dist = 1)
#' plot(rnet_grouped, col = rnet_grouped$group, lwd =
#' rnet_grouped$length / mean(rnet_grouped$length) * 3)
#' routes_fast$group <- rep(1:3, length.out = nrow(routes_fast))
#' rnet_grouped <- overline(routes_fast, attrib = "length", byvars = "group", buff_dist = 1)
#' plot(rnet_grouped,
#' col = rnet_grouped$group, lwd =
#' rnet_grouped$length / mean(rnet_grouped$length) * 3
#' )
#' # sf methods
#' sl = routes_fast_sf[2:4, ]
#' sl <- routes_fast_sf[2:4, ]
#' overline(sl = sl, attrib = "length", buff_dist = 10)
#' rnet_sf = overline(routes_fast_sf, attrib = "length", buff_dist = 10)
#' rnet_sf <- overline(routes_fast_sf, attrib = "length", buff_dist = 10)
#' plot(rnet_sf$geometry, lwd = rnet_sf$length / mean(rnet_sf$length))
#' }
overline <- function(sl, attrib, fun = sum, na.zero = FALSE, byvars = NA, buff_dist = 0) {
UseMethod("overline")
}
#' @export
overline.sf <- function(sl, attrib, fun = sum, na.zero = FALSE, byvars = NA, buff_dist = 0) {

sl_spatial <- sp::SpatialLinesDataFrame(sl = sf::as_Spatial(sl$geometry), data = sf::st_set_geometry(sl, NULL), match.ID = FALSE)
rnet_sp <- overline(sl_spatial, attrib, fun = fun, na.zero = na.zero, byvars = byvars, buff_dist = buff_dist)
sf::st_as_sf(rnet_sp)
@@ -149,94 +149,111 @@ overline.sf <- function(sl, attrib, fun = sum, na.zero = FALSE, byvars = NA, buf
# aggs
}
#' @export
overline.Spatial <- function(sl, attrib, fun = sum, na.zero = FALSE, byvars = NA, buff_dist = 0){

overline.Spatial <- function(sl, attrib, fun = sum, na.zero = FALSE, byvars = NA, buff_dist = 0) {
fun <- c(fun)
if (length(fun) < length(attrib)) {
fun <- rep(c(fun),length.out=length(attrib))
fun <- rep(c(fun), length.out = length(attrib))
}

sl_sp <- as(sl, "SpatialLines")

if(is.na(byvars[1]) == TRUE) {
if (is.na(byvars[1]) == TRUE) {
## get the line sections that make the network
slu <- gsection(sl, buff_dist = buff_dist)
## overlay network with routes
overs = sp::over(slu, sl_sp, returnList = TRUE)
overs <- sp::over(slu, sl_sp, returnList = TRUE)
## overlay is true if end points overlay, so filter them out:
overs = lapply(1:length(overs), function(islu){
Filter(function(isl){
islines(sl_sp[isl,], slu[islu,])
overs <- lapply(1:length(overs), function(islu) {
Filter(function(isl) {
islines(sl_sp[isl, ], slu[islu, ])
}, overs[[islu]])
})
## now aggregate the required attribibute using fun():
#aggs = sapply(overs, function(os){fun(sl[[attrib]][os])})
# aggs = sapply(overs, function(os){fun(sl[[attrib]][os])})
aggs <- setNames(
as.data.frame(
lapply(1:length(attrib),
function(y, overs, attribs, aggfuns){
sapply(overs, function(os,attrib,fun2){
fun2(sl[[attrib]][os])},
attrib=attribs[y],
fun2=aggfuns[[y]])
},
overs,
attrib,
fun)),
attrib)
lapply(
1:length(attrib),
function(y, overs, attribs, aggfuns) {
sapply(overs, function(os, attrib, fun2) {
fun2(sl[[attrib]][os])
},
attrib = attribs[y],
fun2 = aggfuns[[y]]
)
},
overs,
attrib,
fun
)
),
attrib
)

## make a sl with the named attribibute:
sl = sp::SpatialLinesDataFrame(slu, aggs)
#names(sl) = attrib
sl <- sp::SpatialLinesDataFrame(slu, aggs)
# names(sl) = attrib
} else {

splitlines <- lapply(
split(sl, sl@data[,byvars]),
function(x, attrib, gvar){
groupingcat <- unname(unlist(unique(x@data[,gvar])))
sl_spg = as(x, "SpatialLines")
slu = gsection(sl, buff_dist = buff_dist)
overs = sp::over(slu, sl_spg, returnList = TRUE)
overs = lapply(1:length(overs), function(islu) {
Filter(function(isl){islines(sl_spg[isl,],slu[islu,])}, overs[[islu]])
split(sl, sl@data[, byvars]),
function(x, attrib, gvar) {
groupingcat <- unname(unlist(unique(x@data[, gvar])))
sl_spg <- as(x, "SpatialLines")
slu <- gsection(sl, buff_dist = buff_dist)
overs <- sp::over(slu, sl_spg, returnList = TRUE)
overs <- lapply(1:length(overs), function(islu) {
Filter(function(isl) {
islines(sl_spg[isl, ], slu[islu, ])
}, overs[[islu]])
})
#aggs = sapply(overs, function(os){fun(x[[attrib]][os])})
# aggs = sapply(overs, function(os){fun(x[[attrib]][os])})
aggs <- setNames(
as.data.frame(
lapply(1:length(attrib),
function(y, overs, attribs, aggfuns){
sapply(overs, function(os,attrib,fun2){
fun2(x[[attrib]][os])},
attrib=attribs[y],
fun2=aggfuns[[y]])
},
overs,
attrib,
fun)
),
attrib)
sl = sp::SpatialLinesDataFrame(slu, cbind(aggs,as.data.frame(matrix(groupingcat,nrow=1))))
names(sl) = c(attrib,gvar)
sl <- sp::spChFIDs(sl, paste(paste(groupingcat,collapse='.'),row.names(sl@data),sep='.'))
lapply(
1:length(attrib),
function(y, overs, attribs, aggfuns) {
sapply(overs, function(os, attrib, fun2) {
fun2(x[[attrib]][os])
},
attrib = attribs[y],
fun2 = aggfuns[[y]]
)
},
overs,
attrib,
fun
)
),
attrib
)
sl <- sp::SpatialLinesDataFrame(slu, cbind(aggs, as.data.frame(matrix(groupingcat, nrow = 1))))
names(sl) <- c(attrib, gvar)
sl <- sp::spChFIDs(sl, paste(paste(groupingcat, collapse = "."), row.names(sl@data), sep = "."))
sl
},
attrib,
c(byvars)
)

splitlinesdf <- data.frame(dplyr::bind_rows(lapply(splitlines, function(x){x@data})))
row.names(splitlinesdf) <- unname(unlist(lapply(splitlines, function(x){row.names(x@data)})))
splitlinesdf <- data.frame(dplyr::bind_rows(lapply(splitlines, function(x) {
x@data
})))
row.names(splitlinesdf) <- unname(unlist(lapply(splitlines, function(x) {
row.names(x@data)
})))

sl <- sp::SpatialLinesDataFrame(
sp::SpatialLines(unlist(lapply(splitlines, function(x){x@lines}), recursive = FALSE),
proj4string = splitlines[[1]]@proj4string),
sp::SpatialLines(unlist(lapply(splitlines, function(x) {
x@lines
}), recursive = FALSE),
proj4string = splitlines[[1]]@proj4string
),
splitlinesdf
)

}

## remove lines with attribute values of zero
if(na.zero == TRUE){
if (na.zero == TRUE) {
sl <- sl[sl[[attrib]] > 0, ]
}

@@ -288,31 +305,34 @@ onewaygeo <- function(x, attrib) {
onewaygeo.sf <- function(x, attrib) {
geq <- sf::st_equals(x, x, sparse = FALSE) | sf::st_equals_exact(x, x, sparse = FALSE, par = 0.0)
sel1 <- !duplicated(geq) # repeated rows
x$matching_rows = apply(geq, 1, function(x) paste0(formatC(which(x), width = 4, format = "d", flag = 0), collapse = "-"))
x$matching_rows <- apply(geq, 1, function(x) paste0(formatC(which(x), width = 4, format = "d", flag = 0), collapse = "-"))

singlelines <- stats::aggregate(x[attrib], list(x$matching_rows), FUN = sum)

return(singlelines)
}
#' @export
onewaygeo.Spatial <- function(x, attrib){
onewaygeo.Spatial <- function(x, attrib) {
geq <- rgeos::gEquals(x, x, byid = TRUE) | rgeos::gEqualsExact(x, x, byid = TRUE)
sel1 <- !duplicated(geq) # repeated rows
singlelines <- x[sel1,]
singlelines <- x[sel1, ]
non_numeric_cols <- which(!sapply(x@data, is.numeric))
keeper_cols <- sort(unique(c(non_numeric_cols, attrib)))

singlelines@data[, attrib] <- (matrix(
unlist(
lapply(
apply(geq, 1, function(x){
apply(geq, 1, function(x) {
which(x == TRUE)
}),
function(y, x) {
colSums(x[y, attrib]@data)
}, x)),
}, x
)
),
nrow = nrow(x),
byrow=TRUE))[sel1,]
byrow = TRUE
))[sel1, ]

singlelines@data <- singlelines@data[keeper_cols]

32 R/quadrants.R
View file
@@ -11,28 +11,28 @@
#' @export
#' @examples
#' data(zones)
#' sp_obj = zones
#' (quads = quadrant(sp_obj))
#' sp_obj <- zones
#' (quads <- quadrant(sp_obj))
#' plot(sp_obj, col = factor(quads))
#' points(rgeos::gCentroid(sp_obj), col = "white")
#' # edge cases (e.g. when using rasters) lead to NAs
#' sp_obj = raster::rasterToPolygons(raster::raster(ncol = 3, nrow = 3))
#' (quads = quadrant(sp_obj))
#' sp_obj <- raster::rasterToPolygons(raster::raster(ncol = 3, nrow = 3))
#' (quads <- quadrant(sp_obj))
#' plot(sp_obj, col = factor(quads))
quadrant <- function(sp_obj, number_out = FALSE) {
cent = rgeos::gCentroid(sp_obj)
cents = rgeos::gCentroid(sp_obj, byid = TRUE)
in_quadrant = rep(NA, length(sp_obj))
if(number_out) {
in_quadrant[cents@coords[,1] > cent@coords[,1] & cents@coords[,2] > cent@coords[,2]] = 1
in_quadrant[cents@coords[,1] > cent@coords[,1] & cents@coords[,2] < cent@coords[,2]] = 2
in_quadrant[cents@coords[,1] < cent@coords[,1] & cents@coords[,2] > cent@coords[,2]] = 3
in_quadrant[cents@coords[,1] < cent@coords[,1] & cents@coords[,2] < cent@coords[,2]] = 4
cent <- rgeos::gCentroid(sp_obj)
cents <- rgeos::gCentroid(sp_obj, byid = TRUE)
in_quadrant <- rep(NA, length(sp_obj))
if (number_out) {
in_quadrant[cents@coords[, 1] > cent@coords[, 1] & cents@coords[, 2] > cent@coords[, 2]] <- 1
in_quadrant[cents@coords[, 1] > cent@coords[, 1] & cents@coords[, 2] < cent@coords[, 2]] <- 2
in_quadrant[cents@coords[, 1] < cent@coords[, 1] & cents@coords[, 2] > cent@coords[, 2]] <- 3
in_quadrant[cents@coords[, 1] < cent@coords[, 1] & cents@coords[, 2] < cent@coords[, 2]] <- 4
} else {
in_quadrant[cents@coords[,1] > cent@coords[,1] & cents@coords[,2] > cent@coords[,2]] = "NE"
in_quadrant[cents@coords[,1] > cent@coords[,1] & cents@coords[,2] < cent@coords[,2]] = "SE"
in_quadrant[cents@coords[,1] < cent@coords[,1] & cents@coords[,2] > cent@coords[,2]] = "SW"
in_quadrant[cents@coords[,1] < cent@coords[,1] & cents@coords[,2] < cent@coords[,2]] = "NW"
in_quadrant[cents@coords[, 1] > cent@coords[, 1] & cents@coords[, 2] > cent@coords[, 2]] <- "NE"
in_quadrant[cents@coords[, 1] > cent@coords[, 1] & cents@coords[, 2] < cent@coords[, 2]] <- "SE"
in_quadrant[cents@coords[, 1] < cent@coords[, 1] & cents@coords[, 2] > cent@coords[, 2]] <- "SW"
in_quadrant[cents@coords[, 1] < cent@coords[, 1] & cents@coords[, 2] < cent@coords[, 2]] <- "NW"
}
in_quadrant
}
16 R/radiate.R
View file
@@ -29,19 +29,19 @@
#' sum(cents$population) # the total inter-zonal flow
#' plot(flowlines_radiation, lwd = flowlines_radiation$flow / 100)
#' points(cents, cex = cents$population / 100)
od_radiation <- function(p, pop_var = "population", proportion = 1){
od_radiation <- function(p, pop_var = "population", proportion = 1) {
l <- points2flow(p)
l$flow <- NA
for(i in 1:nrow(p)){
for(j in 1:nrow(p)){
if(i == j) next()
for (i in 1:nrow(p)) {
for (j in 1:nrow(p)) {
if (i == j) next()
m <- p[[pop_var]][i]
n <- p[[pop_var]][j]
sel_flow = which(l$O == p@data[i,1] & l$D == p@data[j,1])
sel_flow <- which(l$O == p@data[i, 1] & l$D == p@data[j, 1])
# create circle the radius of which is the distance between i and j centered on i
radius = gprojected(shp = l[sel_flow,], fun = rgeos::gLength)
s_circle <- buff_geo(shp = p[i,], width = radius)
ps <- p[-c(i, j),][s_circle,]
radius <- gprojected(shp = l[sel_flow, ], fun = rgeos::gLength)
s_circle <- buff_geo(shp = p[i, ], width = radius)
ps <- p[-c(i, j), ][s_circle, ]
s <- sum(ps[[pop_var]])
l$flow[sel_flow] <-
p[[pop_var]][i] * proportion * ((m * n) / ((m + s) * (m + n + s)))
35 R/route-transport-api.R
View file
@@ -37,57 +37,60 @@
#' @export
#' @seealso line2route
#' @examples
#'
#'
#' \dontrun{
#' # Plan the 'public' route from Hereford to Leeds
#' rqh <- route_transportapi_public(from = "Hereford", to = "Leeds")
#' plot(rq_hfd)
#' }
#'
# Aim plan public transport routes with transportAPI

#'
#' # Aim plan public transport routes with transportAPI
route_transportapi_public <- function(from, to, silent = FALSE,
region = 'southeast', modes = NA, not_modes = NA){
region = "southeast", modes = NA, not_modes = NA) {

# Convert sp object to lat/lon vector
if(class(from) == "SpatialPoints" | class(from) == "SpatialPointsDataFrame" )
if (class(from) == "SpatialPoints" | class(from) == "SpatialPointsDataFrame") {
from <- coordinates(from)
if(class(to) == "SpatialPoints" | class(to) == "SpatialPointsDataFrame" )
}
if (class(to) == "SpatialPoints" | class(to) == "SpatialPointsDataFrame") {
to <- coordinates(to)
}

# Convert character strings to lon/lat if needs be
if(is.character(from))
if (is.character(from)) {
from <- geo_code(from)
if(is.character(to))
}
if (is.character(to)) {
to <- geo_code(to)
}

orig <- paste0(from, collapse = ",")
dest <- paste0(to, collapse = ",")

api_base = "http://fcc.transportapi.com"
api_base <- "http://fcc.transportapi.com"
ft_string <- paste0("/from/lonlat:", orig, "/to/lonlat:", dest)

queryattrs <- list(region = region)
if (is.na(modes) == FALSE) {
queryattrs[['modes']] = paste0(modes, collapse = "-")
queryattrs[["modes"]] <- paste0(modes, collapse = "-")
} else {
if (is.na(not_modes) == FALSE) {
queryattrs[['not_modes']] = paste0(not_modes, collapse = "-")
queryattrs[["not_modes"]] <- paste0(not_modes, collapse = "-")
}
}

httrreq <- httr::GET(
url = api_base,
path = paste0("/v3/uk/public/journey",ft_string, ".json"),
path = paste0("/v3/uk/public/journey", ft_string, ".json"),
query = queryattrs
)

if (silent == FALSE) {
print(paste0("The request sent to transportapi was: ", httrreq$request$url))
}

if (grepl('application/json',httrreq$headers$`content-type`) == FALSE &
grepl('js',httrreq$headers$`content-type`) == FALSE) {
if (grepl("application/json", httrreq$headers$`content-type`) == FALSE &
grepl("js", httrreq$headers$`content-type`) == FALSE) {
stop("Error: Transportapi did not return a valid result")
}

@@ -106,4 +109,4 @@ route_transportapi_public <- function(from, to, silent = FALSE,

# for the future: add summary data on the route
route
}
}
26 R/route.R
View file
@@ -7,23 +7,24 @@
#' @inheritParams od_coords
#' @inheritParams line2route
#' @export
#' @examples \dontrun{
#' from = "Leek, UK"
#' to = "Hereford, UK"
#' route_leek_to_hereford = route(from, to)
#' @examples
#' \dontrun{
#' from <- "Leek, UK"
#' to <- "Hereford, UK"
#' route_leek_to_hereford <- route(from, to)
#' route(cents_sf[1:3, ], cents_sf[2:4, ]) # sf points
#' route(flowlines_sf[2:4, ]) # lines
#' }
route <- function(from = NULL, to = NULL, l = NULL,
route_fun = route_cyclestreet,
n_print = 10, list_output = FALSE, ...){
route_fun = route_cyclestreet,
n_print = 10, list_output = FALSE, ...) {

# generate od coordinates
FUN <- match.fun(route_fun)
ldf <- od_coords(from, to, l) %>%
dplyr::as_data_frame()

error_fun <- function(e){
error_fun <- function(e) {
warning(paste("Fail for line number", i))
e
}
@@ -40,26 +41,25 @@ route <- function(from = NULL, to = NULL, l = NULL,
rdf[1, ] <- rc[[1]]@data[1, ]
rg[1] <- sf::st_as_sfc(rc[[1]])

if(nrow(ldf) > 1) {
for(i in 2:nrow(ldf)){
if (nrow(ldf) > 1) {
for (i in 2:nrow(ldf)) {
rc[[i]] <- tryCatch({
FUN(from = c(ldf$fx[i], ldf$fy[i]), to = c(ldf$tx[i], ldf$ty[i]), ...)
}, error = error_fun)
perc_temp <- i %% round(nrow(ldf) / n_print)
# print % of distances calculated
if(!is.na(perc_temp) & perc_temp == 0) {
message(paste0(round(100 * i/nrow(ldf)), " % out of ", nrow(ldf), " distances calculated"))
if (!is.na(perc_temp) & perc_temp == 0) {
message(paste0(round(100 * i / nrow(ldf)), " % out of ", nrow(ldf), " distances calculated"))
}

rdf[i, ] <- rc[[i]]@data[1, ]
rg[i] <- sf::st_as_sf(rc[[i]])$geometry

}
}

r <- sf::st_sf(geometry = rg, rdf)

if(list_output) {
if (list_output) {
r <- rc
}