r - How to remove '$osm_points' that are used plot other spatial objects like '$osm_polygons' and '$multipolygons' within osmdata objects in R

Question

I wish find out how many points of interest (POIs, e.g. shops, universities, sports facilities) fall within 400m of various bus stops in Brisbane. In order to do this, I have been extracting features from Open Street Map using the osmdata package in R. My plan is to count the number of points which fall within a 400m buffer of each bus stop, as well as count the number of polygons which intersect with said buffer.

However, I have noticed that some, but not all, of the '$osm_points' within an osmdata object do not seem to represent POIs, but rather the points of more complex spatial objects like '$osm_polygons' and '$osm_multipolygons'. This can be seen in the below example of shops within the Brisbane CBD:

## Import libraries
library(osmdata)
library(ggmap)




# Define functions

## Function to make bounding box
makebb <- function(x) {
  x <- matrix(x,nrow=2,ncol=2)
  rownames(x) <- c("x","y")
  colnames(x) <- c("min","max")
  return(x)
}

## Function to extract all values within a specific key
osm_key <- function(x,bb) {
  x <- bb %>% opq() %>% add_osm_feature(x)
  x <- osmdata_sf(x)
  return(x)
}

## Function to extract specific values within a key
osm_keyvalue <- function(x,y,bb) {
  x <- bb %>% opq() %>% add_osm_feature(x,y)
  x <- osmdata_sf(x)
  return(x)
}

## Function to plot OSM objects
osm_plot <- function(x,bb) {
  defined_map <- get_map(bb, maptype = "toner-background")
  if (is.null(x$osm_multipolygons)) {
    ggmap(defined_map) +
      geom_sf(data = x$osm_points, inherit.aes = FALSE, colour = "#238443",
              fill = "#004529", alpha = 1, size = .7, shape = 21) +
      geom_sf(data = x$osm_polygons, inherit.aes = FALSE, colour = "#238443",
              fill = "#004529", alpha = .5, size = .3, shape = 21) +
      labs(x = "", y = "")
  } else {
    ggmap(defined_map) +
      geom_sf(data = x$osm_points, inherit.aes = FALSE, colour = "#238443",
              fill = "#004529", alpha = 1, size = .7, shape = 21) +
      geom_sf(data = x$osm_polygons, inherit.aes = FALSE, colour = "#238443",
              fill = "#004529", alpha = .5, size = .3, shape = 21) +
      geom_sf(data = x$osm_multipolygons, inherit.aes = FALSE, colour = "#238443",
              fill = "#004529", alpha = .5, size = .3, shape = 21) +
      labs(x = "", y = "")
  }
  
}




# Plot shops in CBD
CBD <- makebb(c(153.02,-27.475,153.03,-27.465))
shop <- osm_key("shop",CBD)
osm_plot(shop,CBD)

As you can see, some points seem to represent independent shops, which I would like to be included in the POI count. However, others seem to simply be the corners of polygons representing larger department stores, which I would like to exclude.

A more extreme example is that of the University of Queensland (UoQ) campus:

# Plot University of Queensland campus
UoQ <- makebb(c(153.00,-27.51,153.03,-27.49))
university <- osm_keyvalue("amenity","university",UoQ)
osm_plot(university,UoQ)

Here, the points clearly don’t represent separate universities, but rather corners of the UoQ campus multipolygon. A bus stop on the edge of campus may therefore incorrectly register several universities within a 400m radius, when there is only one.

Is there any way to remove '$osm_points' which are used to plot other spatial objects like '$osm_polygons' and '$multipolygons', so that only relevant POIs remain within the '$osm_points' section of the osmdata object?

标签： ropenstreetmapintersectionoverlapspatial-query

Answer 1

Probably the best approach here would be creating a new issue in the github issue page of osmdata asking the authors to define a new function to perform the opposite operations of osm_points (or maybe just ask for a suggestion). Meanwhile, I think you can use the following approach:

# packages
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(mapview)

# Define functions

## Function to make bounding box
makebb <- function(x) {
  x <- matrix(x, nrow = 2,ncol = 2)
  rownames(x) <- c("x","y")
  colnames(x) <- c("min","max")
  return(x)
}

## Function to extract all values within a specific key
osm_key <- function(x, bb) {
  x <- bb %>% opq() %>% add_osm_feature(x)
  x <- osmdata_sf(x)
  return(x)
}

## Function to extract specific values within a key
osm_keyvalue <- function(x, y, bb) {
  x <- bb %>% opq() %>% add_osm_feature(x,y)
  x <- osmdata_sf(x)
  return(x)
}

# Plot shops in CBD
CBD <- makebb(c(153.02,-27.475,153.03,-27.465))
shop <- osm_key("shop",CBD)

mapview(shop$osm_points, cex = 2) + shop$osm_polygons

# Get all the points that define the polygons
shop_points_defining_polygons <- osm_points(shop, rownames(shop$osm_polygons))
# Exclude all points that define the polygons
shop_proper_points <- dplyr::anti_join(
  shop$osm_points, 
  sf::st_drop_geometry(shop_points_defining_polygons), 
  by = "osm_id"
)

mapview(shop_proper_points, cex = 2) + shop$osm_polygons

^{Created on 2020-08-23 by the reprex package (v0.3.0)}