###########################################################
# Corrections
library(sf)
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(lwgeom)
## Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.2, PROJ 4.9.3
source("~/Dropbox/BGU/map_projections.R")
setwd("~/Dropbox/Presentations/p_2019_02_LMS_Intro_to_Spatial_R/data_beer_sheva/")
# Fix ITM definition
build = st_read("buildings.geojson", stringsAsFactors = FALSE)
## Reading layer `buildings' from data source `/home/michael/Dropbox/Presentations/p_2019_02_LMS_Intro_to_Spatial_R/data_beer_sheva/buildings.geojson' using driver `GeoJSON'
## replacing null geometries with empty geometries
## Simple feature collection with 23822 features and 8 fields (with 5 geometries empty)
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 34.74723 ymin: 31.20339 xmax: 34.82831 ymax: 31.28896
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
build = st_transform(build, itm3)
build = st_set_crs(build, itm7)
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform
## for that
# Fix topology errors
valid = st_is_valid(build)
all(valid)
## [1] FALSE
build[!valid, ]
## Simple feature collection with 10 features and 8 fields (with 2 geometries empty)
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 177796.1 ymin: 571001.2 xmax: 181538.2 ymax: 577108.7
## epsg (SRID): NA
## proj4string: +proj=tmerc +lat_0=31.7343936111111 +lon_0=35.2045169444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +towgs84=-24.0024,-17.1032,-17.8444,-0.33077,-1.85269,1.66969,5.4248 +units=m +no_defs
## Entrances Not_reside Used floors appartment Name Shape_Leng
## 13891 1 0 מגורים 2 1 <NA> 0.0020807768
## 22003 1 0 מגורים 2 1 <NA> 0.0006821361
## 22005 1 0 מגורים 2 1 <NA> 0.0005988423
## 22104 1 0 מגורים 2 9 <NA> 0.0020103155
## NA NA NA <NA> NA NA <NA> NA
## NA.1 NA NA <NA> NA NA <NA> NA
## 22859 1 0 מגורים 2 24 <NA> 0.0015666372
## 22863 1 0 מגורים 2 24 <NA> 0.0015666372
## 23698 1 0 מגורים 2 7 <NA> 0.0008911881
## 23730 1 0 מגורים 2 14 <NA> 0.0011510526
## Shape_Area geometry
## 13891 3.311376e-08 POLYGON ((178905.1 573177.9...
## 22003 1.854375e-08 POLYGON ((177807.1 571067.7...
## 22005 1.615212e-08 POLYGON ((177805.6 571002.8...
## 22104 8.828333e-08 POLYGON ((181347.7 577108.7...
## NA NA POLYGON EMPTY
## NA.1 NA POLYGON EMPTY
## 22859 4.935403e-08 POLYGON ((179227.5 575104, ...
## 22863 4.935403e-08 POLYGON ((179227.5 575104, ...
## 23698 2.786378e-08 POLYGON ((181533.9 576769.7...
## 23730 4.478178e-08 POLYGON ((181114.1 576604.2...
build = st_make_valid(build)
# Filter empty
empty = st_is_empty(build)
build = build[!empty, ]
# Cast to polygon
build = st_cast(build, "POLYGON")
## Warning in st_cast.MULTIPOLYGON(X[[i]], ...): polygon from first part only
## Warning in st_cast.GEOMETRYCOLLECTION(X[[i]], ...): only first part of
## geometrycollection is retained
# Write
st_write(build, "buildings_corrected.shp", layer_options = "ENCODING=UTF-8", delete_dsn = TRUE)
## Deleting source `buildings_corrected.shp' using driver `ESRI Shapefile'
## Writing layer `buildings_corrected' to data source `buildings_corrected.shp' using driver `ESRI Shapefile'
## options: ENCODING=UTF-8
## features: 23817
## fields: 8
## geometry type: Polygon
###########################################################
# Create apartment layer
library(sf)
library(nngeo)
setwd("~/Dropbox/Presentations/p_2019_02_LMS_Intro_to_Spatial_R/data_beer_sheva/")
# Read layers
build = st_read("buildings_corrected.shp", stringsAsFactors = FALSE)
## Reading layer `buildings_corrected' from data source `/home/michael/Dropbox/Presentations/p_2019_02_LMS_Intro_to_Spatial_R/data_beer_sheva/buildings_corrected.shp' using driver `ESRI Shapefile'
## Simple feature collection with 23817 features and 8 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 175901.6 ymin: 568047.4 xmax: 183624.9 ymax: 577537.9
## epsg (SRID): NA
## proj4string: +proj=tmerc +lat_0=31.7343936111111 +lon_0=35.2045169444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +units=m +no_defs
stat = st_read("statisticalareas_demography2016_WGS84.gpkg", stringsAsFactors = FALSE)
## Reading layer `statisticalareas_demography2016_WGS84' from data source `/home/michael/Dropbox/Presentations/p_2019_02_LMS_Intro_to_Spatial_R/data_beer_sheva/statisticalareas_demography2016_WGS84.gpkg' using driver `GPKG'
## Simple feature collection with 3131 features and 24 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 34.26746 ymin: 29.4905 xmax: 35.89538 ymax: 33.33517
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
# Set building ID
build$b_id = 1:nrow(build)
# Subset
build = build[, c("b_id", "appartment")]
# AOI
ctr = st_centroid(build)
## Warning in st_centroid.sf(build): st_centroid assumes attributes are
## constant over geometries of x
aoi = st_buffer(st_centroid(st_union(ctr)), 250) + c(1500, 3000)
aoi = st_set_crs(aoi, st_crs(build))
# Plot
plot(aoi, border = NA)
plot(st_geometry(build), add = TRUE)
plot(st_geometry(ctr), add = TRUE, col = "red")
# Duplicate
ids = rep(1:nrow(ctr), ctr$appartment)
apartments = ctr[ids, ]
apartments$a_id = 1:nrow(apartments)
#######################################
# Random sampling
set.seed(2)
ids = sample(1:nrow(apartments), 5000)
sample1 = apartments[ids, ]
# Plot
plot(aoi, border = NA)
plot(st_geometry(build), add = TRUE)
plot(st_geometry(sample1), add = TRUE, col = "red")
# How many per building?
sample1$w = 1
build1 = build[sample1, ]
build1 = st_join(build1, sample1[, "w"])
build1 = aggregate(build1[, "w"], st_set_geometry(build1[, "b_id"], NULL), sum)
# Plot
plot(aoi, border = NA)
box()
plot(st_geometry(build), add = TRUE)
plot(st_geometry(build1), add = TRUE, col = "yellow")
text(st_coordinates(st_centroid(build1)), as.character(build1$w), add = TRUE, pos = 3)
## Warning in st_centroid.sf(build1): st_centroid assumes attributes are
## constant over geometries of x
## Warning in text.default(st_coordinates(st_centroid(build1)),
## as.character(build1$w), : "add" is not a graphical parameter
# Filter one per building
d = duplicated(sample1$b_id)
sample2 = sample1[!d, ]
# How many per building?
sample2$w = 1
build2 = build[sample2, ]
build2 = st_join(build2, sample2[, "w"])
build2 = aggregate(build2[, "w"], st_set_geometry(build2[, "b_id"], NULL), sum)
# Plot
plot(aoi, border = NA)
box()
plot(st_geometry(build), add = TRUE)
plot(st_geometry(build2), add = TRUE, col = "yellow")
text(st_coordinates(st_centroid(build2)), as.character(build2$w), add = TRUE, pos = 3)
## Warning in st_centroid.sf(build2): st_centroid assumes attributes are
## constant over geometries of x
## Warning in text.default(st_coordinates(st_centroid(build2)),
## as.character(build2$w), : "add" is not a graphical parameter
#######################################
# Grid-based sampling
# Calculate grid
grid = st_sample(st_buffer(aoi, 5000), type = "regular", size = 10000)
grid = st_set_crs(grid, st_crs(build))
# Select nearest apartment
sample3 = st_join(st_sf(grid), apartments, join = st_nn, maxdist = 50)
sample3 = sample3[!is.na(sample3$a_id), ]
sample3 = apartments[sample3$a_id, ]
# Calculate lines
lines3 = st_connect(sample3, st_sf(grid), progress = FALSE)
# Plot
plot(aoi, border = NA)
plot(st_geometry(build), add = TRUE, border = "grey")
plot(st_geometry(grid), add = TRUE, col = "red")
plot(st_geometry(sample3), add = TRUE, col = "blue")
plot(st_geometry(lines3), add = TRUE, col = "grey50")
text(st_coordinates(st_centroid(lines3)), as.character(round(st_length(lines3))))
# How many per building?
sample3$w = 1
build3 = build[sample3, ]
build3 = st_join(build3, sample3[, "w"])
build3 = aggregate(build3[, "w"], st_set_geometry(build3[, "b_id"], NULL), sum)
# Plot
plot(aoi, border = NA)
box()
plot(st_geometry(build), add = TRUE)
plot(st_geometry(build3), add = TRUE, col = "yellow")
text(st_coordinates(st_centroid(build3)), as.character(build3$w), add = TRUE, pos = 3)
## Warning in st_centroid.sf(build3): st_centroid assumes attributes are
## constant over geometries of x
## Warning in text.default(st_coordinates(st_centroid(build3)),
## as.character(build3$w), : "add" is not a graphical parameter
