points.R

library(units)
library(sf)
library(cartography)

# library(osmdata)
# 
# # define a bounding box
# q0 <- opq(bbox = c(2.2247, 48.8188, 2.4611, 48.9019)) 
# 
# # extract Paris boundaries
# q1 <- add_osm_feature(opq = q0, key = 'name', value = "Paris")
# res1 <- osmdata_sf(q1)
# paris <- st_geometry(res1$osm_multipolygons[1,])
# 
# # extract the Seine river
# q2 <- add_osm_feature(opq = q0, key = 'name', value = "La Seine")
# res2 <- osmdata_sf(q2)
# seine1 <- st_geometry(res2$osm_multilines)
# q2b <- add_osm_feature(opq = q0, key = 'name', 
#                        value = "La Seine - Bras de la Monnaie")
# res2b <- osmdata_sf(q2b)
# seine2 <- st_geometry(res2b$osm_lines)
# 
# # extract Parks and Cemetaries
# q3 <- add_osm_feature(opq = q0, key = 'leisure', value = "park")
# res3 <- osmdata_sf(q3)
# parc1 <- st_geometry(res3$osm_polygons)
# parc2 <- st_geometry(res3$osm_multipolygons)
# q4 <- add_osm_feature(opq = q0, key = 'landuse', value = "cemetery")
# res4 <- osmdata_sf(q4)
# parc3 <- st_geometry(res4$osm_polygons)
# 
# # extract Quartiers
# q5 <- add_osm_feature(opq = q0, key = 'admin_level', value = "10")
# res5 <- osmdata_sf(q5)
# quartier <- res5$osm_multipolygons
# 
# # extract Bars & Pubs
# q6 <- add_osm_feature(opq = q0, key = 'amenity', value = "bar")
# res6 <- osmdata_sf(q6)
# bar <- res6$osm_points
# q7 <- add_osm_feature(opq = q0, key = 'amenity', value = "pub")
# res7 <- osmdata_sf(q7)
# pub <- res7$osm_points
# 
# 
# 
# 
# # use Lambert 93 projection (the french cartographic projection) for all layers
# parc1 <- st_transform(parc1, 2154)
# parc2 <- st_transform(parc2, 2154)
# parc3 <- st_transform(parc3, 2154)
# paris <- st_transform(paris, 2154)
# seine1 <- st_transform(seine1, 2154)
# seine2 <- st_transform(seine2, 2154)
# quartier <- st_transform(quartier, 2154)
# bar <- st_transform(bar, 2154)
# pub <- st_transform(pub, 2154)
# 
# 
# # make layers pretty
# ## Parcs and cemetaries are merged into a single layer, we only keep objects 
# ## greater than 1 ha
# parc <- do.call(c, list(parc1, parc2, parc3))
# parc <- st_union(x = st_buffer(parc,0), by_feature = F)
# parc <- st_cast(parc, "POLYGON")
# parc <- parc[st_area(parc)>=set_units(10000, "m^2")]
# parc <- st_intersection(x = parc, y = paris)
# 
# ## We only keep the part of the river within Paris boundaries
# seine <- st_intersection(x = seine1, y = paris)
# seine <- c(st_cast(seine[1])[2:5], seine[2])
# seine <-c(seine, seine2)
# 
# ## We only keep the bars and pubs within Paris boundaries
# bar <- bar[!is.na(bar$name),]
# pub <- pub[!is.na(pub$name),]
# bars <- rbind(bar[,"name"], pub[,"name"])
# bars <- st_intersection(x = bars, y = paris)
# 
# ## We only keep the Paris quartiers
# quartier <- quartier[substr(quartier$ref.TRIRIS,1,2)==75,]
# 
# 
# save(list= c("paris", "quartier", "seine", "parc", "bars"), 
#      file = "paname.RData", compress = "xz")

load("paname.RData")

grid <- st_make_grid(x = paris, cellsize = 500, square = F)

inter <- st_intersects(x = grid, y = bars)
grid <- st_sf(n = sapply(inter, length), grid)


bars$name

b <- grep("bleu", tolower(bars$name))
b <- c(b, grep("blue", tolower(bars$name)))



r <- grep("red", tolower(bars$name))
r <- c(r, grep("rouge", tolower(bars$name)))
r <- r[-2]

getFigDim(paris, width = 800, mar = c(0,0,1.2,0), res = 100)

png("bars.png", width = 800, height = 450, res = 100)
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier2), col = NA,lty = 2, lwd = 0.2, add=T)
plot(st_geometry(bars), add=T, col = "#0000ff", pch = 20, cex = 0.2)
plot(paris, add=T, lwd = 0.7)
legend(x = "right", legend = c("Bars & Pubs"), 
       col = c("#0000ff"), 
       pch = 20, pt.cex = c(0.2), cex = 0.7, bty = 'n')
layoutLayer(title = "Bars and Pubs in Paris", scale = 1,
            north = T,
            tabtitle = TRUE, frame = FALSE,
            author = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            sources = "T. Giraud, 2020") 
dev.off()


png("bars2.png", width = 800, height = 450, res = 100)
library(cartography)
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "#FBEDDA")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
choroLayer(grid[grid$n>0,], var = "n", border = NA, 
           breaks = c(1,5,10,25,52), 
           add = T, col = hcl.colors(n= 5, "Oslo")[2:5], 
           legend.pos = "right", legend.title.txt = "Number of bars & pubs", legend.horiz = T)
plot(paris, add=T, lwd = 0.7)
layoutLayer(title = "Bars and Pubs in Paris", scale = 1,
            north = T,
            tabtitle = TRUE, frame = FALSE,
            author = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            sources = "T. Giraud, 2020") 
dev.off()





png("bars3.png", width = 800, height = 450, res = 100)
library(cartography)
par(mar = c(0,0,1.2,0))
plot(paris, col = "#D9D0C9", border = NA, bg = "lightblue")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
plot(st_geometry(bars[b,]), add=T, col = "#0000ff", pch = 20, cex = 1)
plot(paris, add=T, lwd = 0.7)
layoutLayer(title = 'Bars and Pubs in Paris with "blue" or "bleu" in their name' , 
            scale = 1,
            north = T,
            tabtitle = TRUE, frame = FALSE,
            author = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            sources = "T. Giraud, 2020") 
dev.off()



png("bars4.png", width = 800, height = 450, res = 100)
library(cartography)
## you need to install the B612 font family
## https://b612-font.com/
par(mar = c(0,0,1.2,0), family = "B612")
plot(paris, col = "#D9D0C9", border = NA, bg = "firebrick1")
plot(parc, col = "#CDEBB2", border = NA, add=T)
plot(seine, col = "#AAD3DF", add=T, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=T)
plot(st_geometry(bars[r,]), add=T, col = "firebrick4", pch = 20, cex = 1)
plot(paris, add=T, lwd = 0.7)
layoutLayer(title = 'Bars and Pubs in Paris with "red" or "rouge" in their name' , 
            scale = 1,
            north = T,
            tabtitle = TRUE, frame = FALSE,
            author = "Map data © OpenStreetMap contributors, under CC BY SA.", 
            sources = "T. Giraud, 2020") 
dev.off()