Detecting streets

library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2

setwd("/home/michael/Dropbox/BGU/Bin_Zhou/p_01_detect_streets_in_building_layer_with_voronoi/")

# Read layer
dat = st_read("BS_Buildings_UTM36N.shp")

## Reading layer `BS_Buildings_UTM36N' from data source `/home/michael/Dropbox/BGU/Bin_Zhou/p_01_detect_streets_in_building_layer_with_voronoi/BS_Buildings_UTM36N.shp' using driver `ESRI Shapefile'
## Simple feature collection with 17460 features and 22 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 665859.9 ymin: 3454151 xmax: 675431.8 ymax: 3463308
## epsg (SRID):    32636
## proj4string:    +proj=utm +zone=36 +datum=WGS84 +units=m +no_defs

# Area of interest
aoi = st_as_sfc("POLYGON ((34.7948 31.2604, 34.7948 31.265, 34.8005 31.265, 34.8005 31.2604, 34.7948 31.2604))")
st_crs(aoi) = 4326
aoi = st_transform(aoi, st_crs(dat))
aoi2 = st_buffer(aoi, 100)

# Clip
dat = dat[aoi2, ]

# Plot 1 - Buildings
plot(dat %>% st_geometry)
plot(aoi %>% st_geometry, border = "yellow", add = TRUE)
plot(aoi2 %>% st_geometry, border = "red", add = TRUE)

# Polygons to lines
dat = st_cast(dat, "MULTILINESTRING")

# Sampling points
pnt = dat %>% st_sample(10000)
pnt = st_sf(pnt)
pnt$id = dat$OBJECTID
pnt_u = st_union(pnt)

# Plot 2 - Sampling points
plot(dat %>% st_geometry)
plot(pnt_u, col = "blue", add = TRUE)
plot(aoi %>% st_geometry, border = "yellow", add = TRUE)
plot(aoi2 %>% st_geometry, border = "red", add = TRUE)

# Voronoi
v = pnt_u %>% st_voronoi %>% st_cast %>% st_intersection(aoi) %>% st_sf

# Plot 3 - Voronoi polygons
plot(v)
plot(aoi %>% st_geometry, border = "yellow", add = TRUE)
plot(aoi2 %>% st_geometry, border = "red", add = TRUE)

# Join with building ID
v = st_join(v, pnt)
v = v[!is.na(v$id), ]

# Aggregate
v = aggregate(v, list(ID = v$id), unique)

# To lines
lines = v %>% st_cast("MULTILINESTRING") %>% st_cast("LINESTRING") %>% st_geometry %>% st_intersection(st_buffer(aoi, dist = -1))

## Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
## geometries for which they may not be constant

# Plot 4 - Final result
plot(dat %>% st_geometry)
plot(lines, col = "blue", add = TRUE)
plot(aoi %>% st_geometry, border = "yellow", add = TRUE)
plot(aoi2 %>% st_geometry, border = "red", add = TRUE)