R package for downloading OpenStreetMap data

Use sf functions instead of sp::point.in.polygon() in trim_osmdata

+ improve behavior: now checks the full geometry instead of just the points
to determine if an object is properly contained in the bbox (only for
osmdata_sf objects, osdmata_sc still wrong).
+ remove reproj form dependencies as it's not used anymore

+168 -178
+1 -2
DESCRIPTION
··· 1 1 Package: osmdata 2 2 Title: Import 'OpenStreetMap' Data as Simple Features or Spatial Objects 3 - Version: 0.2.5.065 3 + Version: 0.2.5.066 4 4 Authors@R: c( 5 5 person("Mark", "Padgham", , "mark.padgham@email.com", role = c("aut", "cre")), 6 6 person("Bob", "Rudis", role = "aut"), ··· 35 35 httr2, 36 36 methods, 37 37 Rcpp (>= 0.12.4), 38 - reproj, 39 38 rvest, 40 39 tibble, 41 40 utils,
+1
NAMESPACE
··· 3 3 S3method(bb_poly_to_mat,SpatialPolygonsDataFrame) 4 4 S3method(bb_poly_to_mat,default) 5 5 S3method(bb_poly_to_mat,list) 6 + S3method(bb_poly_to_mat,matrix) 6 7 S3method(bb_poly_to_mat,sf) 7 8 S3method(bb_poly_to_mat,sfc) 8 9 S3method(c,osmdata)
+2
NEWS.md
··· 13 13 - Deprecate `osmdata_sp` (#372) 14 14 - Pre-prend class names `osmdata_sf` and `osmdata` rather than append; thanks to @agila5 (#373) 15 15 - Add `osmadata_data.frame` class to `osmdata_data_frame()` results (#378) 16 + - Reimplement `trim_osmdata()` using `sf` instead of `sp`. Now, it checks the full geometry instead of just the points 17 + to determine if an object is properly contained in the bbox (only for `osmdata_sf` objects, `osdmata_sc` still wrong) (#382). 16 18 17 19 ## Minor changes 18 20
+145 -165
R/trim-osmdata.R
··· 4 4 #' 5 5 #' @param dat An `osmdata` object returned from [osmdata_sf()] or 6 6 #' [osmdata_sc()]. 7 - #' @param bb_poly A matrix representing a bounding polygon obtained with 8 - #' `getbb (..., format_out = "polygon")` (and possibly selected from 9 - #' resultant list where multiple polygons are returned). 10 - #' @param exclude If TRUE, objects are trimmed exclusively, only retaining those 7 + #' @param bb_poly An `sf` or `sfc` object, or matrix representing a bounding 8 + #' polygon. Can be obtained with `getbb (..., format_out = "polygon")` or 9 + #' `getbb (..., format_out = "sf_polygon")` (and possibly 10 + #' selected from resultant list where multiple polygons are returned). 11 + #' @param exclude If `TRUE`, objects are trimmed exclusively, only retaining those 11 12 #' strictly within the bounding polygon; otherwise all objects which partly 12 13 #' extend within the bounding polygon are retained. 13 14 #' ··· 20 21 #' `getbb(..., format_out = "polygon"|"sf_polygon")`. These shapes can be 21 22 #' outdated and thus could cause the trimming operation to not give results 22 23 #' expected based on the current state of the OSM data. 24 + #' @note To reduce the downloaded data from Overpass, you can do the trimming in 25 + #' the server-side using `getbb(..., format_out = "osm_type_id")` 26 + #' (see examples). 23 27 #' 24 28 #' @family transform 25 29 #' ··· 41 45 #' bb_sp <- as (bb_sf, "Spatial") 42 46 #' class (bb_sp) # SpatialPolygonsDataFrame 43 47 #' dat_tr <- trim_osmdata (dat, bb_sp) 48 + #' 49 + #' # Server-side trimming equivalent 50 + #' bb <- getbb ("colchester uk", format_out = "osm_type_id") 51 + #' query <- opq (bb) |> 52 + #' add_osm_feature (key = "highway") 53 + #' dat <- osmdata_sf (query, quiet = FALSE) 44 54 #' } 45 55 #' @export 46 56 trim_osmdata <- function (dat, bb_poly, exclude = TRUE) { ··· 65 75 trim_osmdata.osmdata_sf <- function (dat, bb_poly, exclude = TRUE) { 66 76 67 77 requireNamespace ("sf") 68 - if (!is (bb_poly, "matrix")) { 69 - bb_poly <- bb_poly_to_mat (bb_poly) 78 + if (!inherits (bb_poly, c ("sf", "sfc"), which = FALSE)) { 79 + bb_poly <- bb_poly_to_sf (bb_poly) 70 80 } 71 81 72 - if (nrow (bb_poly) > 1) { 82 + dat <- trim_to_poly_pts (dat, bb_poly, exclude = exclude) 83 + dat <- trim_to_poly (dat, bb_poly = bb_poly, exclude = exclude) 84 + dat <- trim_to_poly_multi (dat, bb_poly = bb_poly, exclude = exclude) 73 85 74 - # "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0": 75 - srcproj <- .lonlat () 76 - # "+proj=merc +a=6378137 +b=6378137": 77 - crs <- .sph_merc () 78 - bb_poly <- reproj::reproj ( 79 - bb_poly, 80 - target = crs, 81 - source = srcproj 82 - ) [, 1:2] 83 - 84 - dat <- trim_to_poly_pts (dat, bb_poly, exclude = exclude) 85 - dat <- trim_to_poly (dat, bb_poly = bb_poly, exclude = exclude) 86 - dat <- trim_to_poly_multi (dat, bb_poly = bb_poly, exclude = exclude) 87 - } else { 88 - message ( 89 - "bb_poly must be a matrix with > 1 row; ", 90 - " data will not be trimmed." 91 - ) 92 - } 93 86 return (dat) 94 87 } 95 88 ··· 124 117 stop ("bb_poly is of unknown class; please use matrix or a spatial class") 125 118 } 126 119 120 + 127 121 more_than_one <- function () { 128 122 129 123 message ("bb_poly has more than one polygon; the first will be selected.") 124 + } 125 + 126 + #' @export 127 + bb_poly_to_mat.matrix <- function (x) { 128 + x 130 129 } 131 130 132 131 #' @export ··· 178 177 return (x) 179 178 } 180 179 181 - trim_to_poly_pts <- function (dat, bb_poly, exclude = TRUE) { 182 180 183 - if (inherits (dat$osm_points, "sf")) { 181 + bb_poly_to_sf <- function (bb_poly) { 184 182 185 - requireNamespace ("sp", quietly = TRUE) ## TODO: remove sp 183 + bb_poly <- bb_poly_to_mat (bb_poly) 186 184 187 - # "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0": 188 - srcproj <- .lonlat () 189 - # "+proj=merc +a=6378137 +b=6378137": 190 - crs <- .sph_merc () 185 + if (nrow (bb_poly) == 2) { # bbox corners 191 186 192 - g <- do.call (rbind, dat$osm_points$geometry) 193 - ## TODO: don't reproject. s2 used by sf works fine with lon lat 194 - # See also bb_poly_* 195 - g <- reproj::reproj (g, target = crs, source = srcproj) 196 - indx <- sp::point.in.polygon ( # TODO: use sf::st_intersects() 197 - g [, 1], g [, 2], 198 - bb_poly [, 1], bb_poly [, 2] 187 + bb_poly <- rbind ( 188 + bb_poly [1, ], 189 + c (bb_poly [1, 1], bb_poly [2, 2]), 190 + bb_poly [2, ], 191 + c (bb_poly [2, 1], bb_poly [1, 2]) 199 192 ) 193 + } 194 + 195 + if (!identical ( 196 + as.numeric (utils::head (bb_poly, 1)), 197 + as.numeric (utils::tail (bb_poly, 1)) 198 + )) { # Non-closed polygon 199 + bb_poly <- rbind (bb_poly, bb_poly [1, ]) 200 + } 201 + 202 + bb_poly <- sf::st_sfc (sf::st_polygon (list (bb_poly)), crs = 4326) 203 + 204 + return (bb_poly) 205 + } 206 + 207 + 208 + #' Remove points outside bb_poly 209 + #' 210 + #' @param dat an object of class `osmdata_sf` 211 + #' @param bb_poly an object of class `sf` or `sfc` 212 + #' @param exclude If TRUE, only retaini points strictly within the bounding 213 + #' polygon and discard points in the vertex or boundaries; otherwise keep the 214 + #' points in the boundaries of the bb_poly. 215 + #' 216 + #' @noRd 217 + trim_to_poly_pts <- function (dat, bb_poly, exclude = TRUE) { 218 + 219 + if (inherits (dat$osm_points, "sf")) { 220 + 200 221 if (exclude) { 201 - indx <- which (indx == 1) 222 + # st_contains_properly assumes planar coordinates. No s2 method? 223 + g <- sf::st_transform (dat$osm_points, crs = .sph_merc ()) 224 + bb_poly <- sf::st_transform (bb_poly, crs = .sph_merc ()) 225 + indx <- sf::st_contains_properly (bb_poly, g, sparse = FALSE) [1, ] 202 226 } else { 203 - indx <- which (indx > 0) 227 + indx <- sf::st_within ( 228 + dat$osm_points, bb_poly, 229 + sparse = FALSE 230 + ) [, 1] 204 231 } 232 + 205 233 dat$osm_points <- dat$osm_points [indx, ] 206 234 } 207 235 208 236 return (dat) 209 237 } 238 + 210 239 211 240 #' get_trim_indx 212 241 #' ··· 223 252 #' @noRd 224 253 get_trim_indx <- function (g, bb, exclude) { 225 254 226 - # "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0": 227 - srcproj <- .lonlat () 228 - # "+proj=merc +a=6378137 +b=6378137": 229 - crs <- .sph_merc () 230 - 231 - indx <- lapply (g, function (i) { 232 - 233 - if (is.list (i)) { # polygons 234 - i <- i [[1]] 235 - } 236 - i <- reproj::reproj (as.matrix (i), target = crs, source = srcproj) 237 - inp <- sp::point.in.polygon ( 238 - i [, 1], i [, 2], 239 - bb [, 1], bb [, 2] 240 - ) 241 - if ((exclude & all (inp > 0)) | 242 - (!exclude & any (inp > 0))) { 243 - return (TRUE) 244 - } else { 245 - return (FALSE) 246 - } 247 - }) 248 - ret <- NULL # multi objects can be empty 249 - if (length (indx) > 0) { 250 - ret <- which (unlist (indx)) 255 + if (exclude) { 256 + # st_contains_properly assumes planar coordinates. No s2 method? 257 + g <- sf::st_transform (g, crs = .sph_merc ()) 258 + bb <- sf::st_transform (bb, crs = .sph_merc ()) 259 + indx <- sf::st_contains_properly (bb, g, sparse = FALSE) [1, ] 260 + } else { 261 + indx <- sf::st_intersects (bb, g, sparse = FALSE) [1, ] 251 262 } 252 - return (ret) 263 + 264 + return (which (indx)) 253 265 } 254 266 255 267 trim_to_poly <- function (dat, bb_poly, exclude = TRUE) { 256 268 257 - if (is (dat$osm_lines, "sf") | is (dat$osm_polygons, "sf")) { 269 + if (inherits (dat$osm_lines, "sf")) { 270 + gnms <- "osm_lines" 271 + } else { 272 + gnms <- character () 273 + } 274 + if (inherits (dat$osm_polygons, "sf")) { 275 + gnms <- c (gnms, "osm_polygons") 276 + } 258 277 259 - gnms <- c ("osm_lines", "osm_polygons") 260 - index <- vapply ( 261 - gnms, function (i) !is.null (dat [[i]]), 262 - logical (1) 263 - ) 264 - gnms <- gnms [index] 265 - for (g in gnms) { 278 + for (g in gnms) { 266 279 267 - if (!is.null (dat [[g]]) & nrow (dat [[g]]) > 0) { 280 + if (nrow (dat [[g]]) > 0) { 268 281 269 - indx <- get_trim_indx (dat [[g]]$geometry, bb_poly, 270 - exclude = exclude 271 - ) 272 - # cl <- class (dat [[g]]$geometry) # TODO: Delete 273 - attrs <- attributes (dat [[g]]) 274 - attrs$row.names <- attrs$row.names [indx] 275 - attrs_g <- attributes (dat [[g]]$geometry) 276 - attrs_g$names <- attrs_g$names [indx] 277 - dat [[g]] <- dat [[g]] [indx, ] # this strips sf class defs 278 - # class (dat [[g]]$geometry) <- cl # TODO: Delete 279 - attributes (dat [[g]]) <- attrs 280 - attributes (dat [[g]]$geometry) <- attrs_g 281 - } 282 + indx <- get_trim_indx ( 283 + g = dat [[g]]$geometry, bb = bb_poly, exclude = exclude 284 + ) 285 + dat [[g]] <- dat [[g]] [indx, ] 282 286 } 283 287 } 284 288 ··· 287 291 288 292 trim_to_poly_multi <- function (dat, bb_poly, exclude = TRUE) { 289 293 290 - if (is (dat$osm_multilines, "sf") | is (dat$osm_multipolygons, "sf")) { 294 + if (inherits (dat$osm_multilines, "sf")) { 295 + gnms <- "osm_multilines" 296 + } else { 297 + gnms <- character () 298 + } 299 + if (inherits (dat$osm_multipolygons, "sf")) { 300 + gnms <- c (gnms, "osm_multipolygons") 301 + } 291 302 292 - gnms <- c ("osm_multilines", "osm_multipolygons") 293 - index <- vapply ( 294 - gnms, function (i) !is.null (dat [[i]]), 295 - logical (1) 296 - ) 297 - gnms <- gnms [index] 298 - for (g in gnms) { 303 + for (g in gnms) { 299 304 300 - if (nrow (dat [[g]]) > 0) { 305 + if (nrow (dat [[g]]) > 0) { 301 306 302 - if (g == "osm_multilines") { 303 - indx <- lapply (dat [[g]]$geometry, function (gi) { 304 - get_trim_indx ( 305 - g = gi, bb = bb_poly, 306 - exclude = exclude 307 - ) 308 - }) 309 - } else { 310 - indx <- lapply (dat [[g]]$geometry, function (gi) { 311 - get_trim_indx ( 312 - g = gi [[1]], bb = bb_poly, 313 - exclude = exclude 314 - ) 315 - }) 316 - } 317 - ilens <- vapply (indx, length, 1L, USE.NAMES = FALSE) 318 - glens <- vapply (dat [[g]]$geometry, function (i) { 319 - length (i [[1]]) 320 - }, 1L, USE.NAMES = FALSE) 321 - if (exclude) { 322 - indx <- which (ilens == glens) 323 - } else { 324 - indx <- which (ilens > 0) 325 - } 307 + # if (g == "osm_multilines") { 308 + # indx <- lapply (dat [[g]]$geometry, function (gi) { 309 + # get_trim_indx ( 310 + # g = gi, bb = bb_poly, 311 + # exclude = exclude 312 + # ) 313 + # }) 314 + # } else { 315 + # indx <- lapply (dat [[g]]$geometry, function (gi) { 316 + # get_trim_indx ( 317 + # g = gi [[1]], bb = bb_poly, 318 + # exclude = exclude 319 + # ) 320 + # }) 321 + # } 322 + indx <- get_trim_indx ( 323 + g = dat [[g]]$geometry, bb = bb_poly, exclude = exclude 324 + ) 325 + # ilens <- vapply (indx, length, 1L, USE.NAMES = FALSE) 326 + # glens <- vapply (dat [[g]]$geometry, function (i) { 327 + # length (i [[1]]) 328 + # }, 1L, USE.NAMES = FALSE) 329 + # if (exclude) { 330 + # indx <- which (ilens == glens) 331 + # } else { 332 + # indx <- which (ilens > 0) 333 + # } 326 334 327 - # cl <- class (dat [[g]]$geometry) # TODO: Delete 328 - attrs <- attributes (dat [[g]]) 329 - attrs$row.names <- attrs$row.names [indx] 330 - attrs_g <- attributes (dat [[g]]$geometry) 331 - attrs_g$names <- attrs_g$names [indx] 332 - dat [[g]] <- dat [[g]] [indx, ] 333 - # class (dat [[g]]$geometry) <- cl # TODO: Delete 334 - attributes (dat [[g]]) <- attrs 335 - attributes (dat [[g]]$geometry) <- attrs_g 336 - } 335 + dat [[g]] <- dat [[g]] [indx, ] 337 336 } 338 337 } 339 338 ··· 347 346 #' @export 348 347 trim_osmdata.osmdata_sc <- function (dat, bb_poly, exclude = TRUE) { 349 348 350 - v <- verts_in_bpoly (dat, bb_poly) 349 + v <- verts_in_bpoly (dat, bb_poly, exclude = exclude) 350 + # TODO: no geometries checked, only vertex 351 351 352 352 if (exclude) { 353 353 ··· 390 390 return (dat) 391 391 } 392 392 393 - verts_in_bpoly <- function (dat, bb_poly) { 393 + verts_in_bpoly <- function (dat, bb_poly, exclude) { 394 394 395 - bb_poly_to_sf <- function (bb_poly) { 396 - 397 - if (nrow (bb_poly) == 2) { 398 - 399 - bb_poly <- rbind ( 400 - bb_poly [1, ], 401 - c (bb_poly [1, 1], bb_poly [2, 2]), 402 - bb_poly [2, ], 403 - c (bb_poly [2, 1], bb_poly [1, 2]) 404 - ) 405 - } 406 - 407 - if (!identical ( 408 - as.numeric (utils::head (bb_poly, 1)), 409 - as.numeric (utils::tail (bb_poly, 1)) 410 - )) { 411 - bb_poly <- rbind (bb_poly, bb_poly [1, ]) 412 - } 413 - 414 - sf::st_sf ( 415 - sf::st_sfc ( 416 - sf::st_polygon (list (bb_poly)), 417 - crs = 4326 418 - ) 419 - ) 420 - } 395 + requireNamespace ("sf") 421 396 bb_poly <- bb_poly_to_sf (bb_poly) 422 397 423 398 vert_to_sf <- function (dat) { ··· 426 401 sf::st_as_sf (v, coords = c ("x_", "y_"), crs = 4326) 427 402 } 428 403 429 - # suppress message about st_intersection assuming planar coordinates, 430 - # because the inaccuracy may be ignored here 431 - suppressMessages (w <- sf::st_within (vert_to_sf (dat), bb_poly)) 404 + if (exclude) { 405 + # st_contains_properly assumes planar coordinates. No s2 method? 406 + g <- sf::st_transform (vert_to_sf (dat), crs = .sph_merc ()) 407 + bb_poly <- sf::st_transform (bb_poly, crs = .sph_merc ()) 408 + w <- sf::st_contains_properly (bb_poly, g, sparse = FALSE) [1, ] 409 + } else { 410 + w <- sf::st_within (vert_to_sf (dat), bb_poly, sparse = FALSE) [, 1] 411 + } 432 412 433 - dat$vertex$vertex_ [which (as.logical (w))] 413 + return (dat$vertex$vertex_ [which (w)]) 434 414 }
+15 -4
man/trim_osmdata.Rd
··· 10 10 \item{dat}{An \code{osmdata} object returned from \code{\link[=osmdata_sf]{osmdata_sf()}} or 11 11 \code{\link[=osmdata_sc]{osmdata_sc()}}.} 12 12 13 - \item{bb_poly}{A matrix representing a bounding polygon obtained with 14 - \code{getbb (..., format_out = "polygon")} (and possibly selected from 15 - resultant list where multiple polygons are returned).} 13 + \item{bb_poly}{An \code{sf} or \code{sfc} object, or matrix representing a bounding 14 + polygon. Can be obtained with \code{getbb (..., format_out = "polygon")} or 15 + \code{getbb (..., format_out = "sf_polygon")} (and possibly 16 + selected from resultant list where multiple polygons are returned).} 16 17 17 - \item{exclude}{If TRUE, objects are trimmed exclusively, only retaining those 18 + \item{exclude}{If \code{TRUE}, objects are trimmed exclusively, only retaining those 18 19 strictly within the bounding polygon; otherwise all objects which partly 19 20 extend within the bounding polygon are retained.} 20 21 } ··· 33 34 \code{getbb(..., format_out = "polygon"|"sf_polygon")}. These shapes can be 34 35 outdated and thus could cause the trimming operation to not give results 35 36 expected based on the current state of the OSM data. 37 + 38 + To reduce the downloaded data from Overpass, you can do the trimming in 39 + the server-side using \code{getbb(..., format_out = "osm_type_id")} 40 + (see examples). 36 41 } 37 42 \examples{ 38 43 \dontrun{ ··· 52 57 bb_sp <- as (bb_sf, "Spatial") 53 58 class (bb_sp) # SpatialPolygonsDataFrame 54 59 dat_tr <- trim_osmdata (dat, bb_sp) 60 + 61 + # Server-side trimming equivalent 62 + bb <- getbb ("colchester uk", format_out = "osm_type_id") 63 + query <- opq (bb) |> 64 + add_osm_feature (key = "highway") 65 + dat <- osmdata_sf (query, quiet = FALSE) 55 66 } 56 67 } 57 68 \seealso{
+4 -7
tests/testthat/test-trim.R
··· 13 13 trim_osmdata (1, bb_poly = bb), 14 14 "unrecognised dat class: numeric" 15 15 ) 16 - expect_silent (x1 <- trim_osmdata (x0, bb_poly = bb)) 16 + expect_silent (x1 <- trim_osmdata (x0, bb_poly = bb, exclude = TRUE)) 17 17 expect_equal (nrow (x1$osm_points), 0) 18 18 expect_equal (nrow (x1$osm_lines), 0) 19 19 expect_equal (nrow (x1$osm_polygons), 0) 20 20 expect_equal (nrow (x1$osm_multilines), 0) 21 21 expect_equal (nrow (x1$osm_multipolygons), 0) 22 22 23 - expect_silent (x1 <- trim_osmdata (x0, 24 - bb_poly = bb, 25 - exclude = FALSE 26 - )) 27 - expect_equal (nrow (x1$osm_points), 2) 23 + expect_silent (x1 <- trim_osmdata (x0, bb_poly = bb, exclude = FALSE)) 24 + expect_equal (nrow (x1$osm_points), 4) 28 25 expect_equal (nrow (x1$osm_lines), 1) 29 26 expect_equal (nrow (x1$osm_polygons), 1) 30 27 expect_equal (nrow (x1$osm_multilines), 0) ··· 81 78 82 79 bb_sf <- sf::st_polygon (list (bb)) |> 83 80 st_sfc () |> 84 - st_sf () 81 + st_sf (crs = 4326) 85 82 expect_silent (x2 <- trim_osmdata (x0, 86 83 bb_poly = bb_sf, 87 84 exclude = FALSE