---
title: "Using OSM Building Data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using OSM Building Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
set.seed(19750604)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

In areas where good OpenStreetMap building surveys are available, this data makes a 
good proxy for a population density when no census is available.
The `osmdata` package lets a user query the OSM database for particular features.
This vignette illustrates how to get OSM building data and use it for a sampling
design.

## Input Polygon

First we need to define a polygon for our study region. Here we construct it from 
a WKT string, but typically this might be read in from a data source such as a 
geopackage or a shapefile, or even drawn on a map using R's
interactove `locator` function.

```{r setup}
library(sidclosepairs)
library(sf)
study_area = "POLYGON (
 (-2.794796 54.02322, -2.801250 54.02686, -2.802311 54.02867, -2.802351 54.02891, -2.802334 54.03084,
  -2.802218 54.03190, -2.801649 54.03199, -2.801142 54.03202, -2.784621 54.03237, -2.780087 54.03196,
  -2.778327 54.03087, -2.777984 54.03048, -2.775405 54.02391, -2.775982 54.02327, -2.794796 54.02322))"

## remove newline markers
study_area = gsub("\n","", study_area)

# convert WKT to spatial polygon data frame
study_area = st_as_sf(data.frame(wkt=study_area), wkt="wkt", crs=4326)
```

## Get the Buildings

Now with our study area defined, we construct a query using `osmdata` package functions and
get back the points in the bounding box of our study area. Taking the centroid of the 
buildings makes sure we reduce polygons down to points. Finally we take only the points
that are within the polygon itself, not just the bounding box. We'll wrap these calls
into a function for re-use and call it to get buildings in our study area.

```{r getem}
library(osmdata)
## or use another overpass server
set_overpass_url("https://overpass-api.de/api/interpreter")
get_buildings <- function(area){
    osmd <- opq(bbox = st_bbox(area))
    osmd <- add_osm_feature(osmd, key = "building")
    osmd <- osmdata_sf(osmd)
    centroids_in_area <- st_intersection(st_centroid(osmd$osm_polygons), area)
    return(centroids_in_area)
}
buildings = get_buildings(study_area)
```

If you want to do this for several study areas you should consider wrapping
the `osmdata` calls into a single function for convenience.

## Inhibitory Close Pair Sampling

Now we'll sample fifty points in total, with ten sets of close pairs, keeping our 
inhibitory points 100 metres apart, and our close pairs in a 20 metre radius.

First we have to convert our data to a local Cartesian coordinate system,
so we use EPSG code 27700, for the Ordnance Survey GB system which is appropriate
for this region in Lancaster.

```{r convo}
crs = 27700
study_gb = st_transform(study_area, crs)
buildings = st_transform(buildings, crs)
```

Now we can sample our 50 points with 10 close pairs using the "census" proposal method
for the inhibitory points, and a 20 metre radius for the close pairs.

```{r samplez, fig.width=6}
pts_buildings = icpSample(50, 10, 100, 20, censusxy=buildings, poly=study_gb, proposal="census")
pts_buildings
```

## Plotting

To show these points we can first plot the building locations in grey, then overlay the
first 40 points in red to show the inhibitory points, and then add the close pairs in blue.

```{r plots, fig.width=6}
plot(study_gb)
plot(st_geometry(buildings), add=TRUE, col="grey", cex=.5)
plot(st_geometry(pts_buildings[1:40,]), add=TRUE, col="red")
plot(st_geometry(pts_buildings[41:50,]), add=TRUE, col="blue")
```

Note that the blue close pair points are not guaranteed to be on a building, but some
other method will need to be used to select the second building for a survey related
to that point. Often the nearest building to the point will be sampled, or the nearest in 
the direction indicated by the two points.