Εντοπισμός της εμπορικής περιοχής σε μια πόλη χρησιμοποιώντας δεδομένα OpenStreetMap. Θα μελετήσουμε το πλήθος, το είδος επιχειρήσεων και τη χωρική κατανομή τους, σε μια υπο-περιοχή της πόλης της Μυτιλήνης, στο Δήμο Μυτιλήνης, στη Λέσβο.
library(remotes)
library(ggplot2)
library(curl)
library(dplyr)
library(sf)
library(terra)
library(devtools)
#devtools::install_github("dimitrisk/goal", quiet=T)
library(goal)
library(osmdata)
# remotes::install_github('ropensci/osmdata')
library(raster)
library(terra)
library(knitr)
library(ggplot2)
opts_chunk$set(cache=TRUE)
options(scipen=999)Διαβάζουμε το πολύγωνο της ακτογραμμής της Λέσβου. Στη συνέχεια ορίζουμε ένα δικό μας τετράγωνο περίγραμμα (Bounding Box) το οποίο θα χρησιμοποιήσουμε για να κατεβάσουμε δεδομένα.
Διαθέσημες πόλεις για ανάλυση (bounding box)
bb_mytilini = c(26.5392, 39.0806, 26.5689, 39.123)
bb_chios = c(26.1092, 38.333, 26.1452, 38.4268)
bb_rodos = c(28.193, 36.4093, 28.2441, 36.4617)
bb_chania = c(23.9639, 35.4894, 24.0611, 35.5322)
bb_patra = c(21.6662, 38.1726, 21.7988, 38.2858)
bb_larisa = c(22.3707, 39.5915, 22.4737, 39.6595)
bb1 = bb_mytilini
bb2 = bb_chiosΜεγεθος raster (EPSG:4326)
Χρησιμοποιώντας ένα Bounding-Box (ΒΒ) κατεβάζουμε όλες τις επιχειρήσεις από το OpenStreetMap. (Points, Polygons, Multipolygons) οι οποίες βρίσκονται εντός του Bounding-Box. Στη συνέχεια συνενώνουμε όλα τα είδη των επιχειρήσεων σε ένα σημειακό σύνολο-δεδομένων.
Simple feature collection with 950 features and 4 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 26.54471 ymin: 39.08545 xmax: 26.56383 ymax: 39.11596
Geodetic CRS: WGS 84
First 10 features:
osm_id name shop geotype geometry
232597560 232597560 Κωτσόβολος electronics frompolygon POINT (26.55529 39.10176)
252369882 252369882 Μy Market supermarket frompolygon POINT (26.55517 39.10033)
255945997 255945997 Κοινωνικό παντοπωλείο supermarket frompolygon POINT (26.55563 39.10689)
258057913 258057913 Μολυβιάτης toys frompolygon POINT (26.55741 39.10794)
258057914 258057914 Μολυβιάτης toys frompolygon POINT (26.55747 39.10811)
258566095 258566095 My Market supermarket frompolygon POINT (26.55724 39.10899)
262504954 262504954 Κινέζικο clothes frompolygon POINT (26.55595 39.10205)
262504958 262504958 ΑΒ Βασιλόπουλος supermarket frompolygon POINT (26.55528 39.10208)
269803274 269803274 Public electronics frompolygon POINT (26.55789 39.1069)
269809517 269809517 Andriana Tours travel_agency frompolygon POINT (26.55849 39.10656)
Μετατρέπουμε το Bounding-Box (ΒΒ) των επιχειρήσεων σε πολύγονο. Υπολογίζουμε την τομή (intersection) μεταξύ ακτογραμμής Λέσβου και Bounding-Box των επιχειρήσεων.
mypol_bb = osm.osmdata_result_2_bbox_pol(shops_all) %>% st_transform("EPSG:4326") # get polygon of this bounding box.
pol = sf::st_intersection( mypol_bb, acto )Δεδομένα μέχρι στιγμής:
Κατασκευάζουμε ένα άδειο raster για χρήση στη συνέχεια.
class : RasterLayer
dimensions : 21, 15, 315 (nrow, ncol, ncell)
resolution : 0.002, 0.002 (x, y)
extent : 26.5392, 26.5692, 39.0806, 39.1226 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
Μετατροπή όλων των δεδομένων σε wgs:2100
Πίνακας συχνοτήτων των επιχειρήσεων στη περιοχή μελέτης.
# A tibble: 92 × 3
shop n freq
* <chr> <int> <dbl>
1 <NA> 341 35.9
2 clothes 91 9.58
3 convenience 49 5.16
4 hairdresser 38 4
5 jewelry 26 2.74
6 shoes 21 2.21
7 bakery 20 2.11
8 kiosk 19 2
9 butcher 16 1.68
10 supermarket 16 1.68
# ℹ 82 more rows
# A tibble: 92 × 3
shop n freq
* <chr> <int> <dbl>
1 <NA> 341 35.9
2 clothes 91 9.58
3 convenience 49 5.16
4 hairdresser 38 4
5 jewelry 26 2.74
6 shoes 21 2.21
7 bakery 20 2.11
8 kiosk 19 2
9 butcher 16 1.68
10 supermarket 16 1.68
11 travel_agency 15 1.58
12 florist 13 1.37
13 confectionery 12 1.26
14 optician 11 1.16
15 alcohol 10 1.05
16 electronics 10 1.05
17 books 8 0.842
18 cosmetics 8 0.842
19 greengrocer 8 0.842
20 coffee 7 0.737
21 computer 7 0.737
22 lottery 7 0.737
23 seafood 7 0.737
24 stationery 7 0.737
25 toys 7 0.737
26 bicycle 6 0.632
27 car_repair 6 0.632
28 mobile_phone 6 0.632
29 yes 6 0.632
30 beauty 5 0.526
# ℹ 62 more rows
Επιλογή μόνο των εμπορικών επιχειρήσεων
emporikes_epixeiriseis = c("clothes","supermarket","bakery","mobile_phone","computer","convenience","hairdresser","jewelry","shoes","florist","electronics","coffee")
shops_ena = shops[shops$shop %in% emporikes_epixeiriseis,]Συχνότητα εμπορικών επιχειρήσεων
Var1 Freq
1 clothes 91
2 convenience 49
3 hairdresser 38
4 jewelry 26
5 shoes 21
6 bakery 20
7 supermarket 16
8 florist 13
9 electronics 10
10 coffee 7
11 computer 7
12 mobile_phone 6
Συνολικό πλήθος εμπορικών επιχειρήσεων
[1] 304
Raster με πλήθος εμπορικών επιχειρήσεων αν κελί
Raster με Ύπαρξη/Απουσία επιχειρήσεων ανά κελί
Μεθοδος για τον υπολογισμό του πινακα αποστάσεων μεταξύ των σημείων. Ειτε για ολα τα σημεια, ειτε για μια κατηγορία σημειων μονο.
myDist = function(inPoints, selecttionColumn="", SelectionValue=""){
if(nchar(selecttionColumn)>1){
myFiltered = inPoints %>% filter(.data[[selecttionColumn]] == SelectionValue)
}else{
myFiltered = inPoints
}
jim = st_distance(myFiltered) %>% units::drop_units() %>% as.matrix()
jim[lower.tri(jim,diag=T)] = NA
return(jim)
}Ποιες ειναι οι αποστάσεις μεταξύ ολων των εμπορικών επιχειρησεων γενικά ?
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.63 205.02 541.01 685.26 1030.70 3388.63 46360
Ποιες ειναι οι αποστάσεις μεταξύ των supermarket ?
myDist(shops_ena, selecttionColumn="shop", SelectionValue="supermarket") %>% as.vector( ) %>% summary() Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2.396 646.353 1076.043 1174.500 1669.772 3002.308 136
Οπτικοποιηση
Πίνακας συχνοτήτων κελιών
1 2 3 4 5 6 8 9 12 19 36 62 65
24 14 3 4 2 1 1 1 1 1 1 1 1
Ιστόγραμμα με το πλήθος εμπορικών επιχειρήσεων ανά κελί για την Μυτιλήνη.
data_hist = as.data.frame(density1)
data_hist$Perioxi = "Mytilini"
p = data_hist %>%
ggplot( aes(x=layer, fill=Perioxi)) +
ggtitle("Ιστόγραμμα συχνοτήτων (πλήθος επιχειρήσεων ανά κελί)")+
geom_histogram( alpha=0.6, position = 'identity') +
xlab("Εμπορικές επιχειρήσεις")+ylab("Συχνότητα")+ xlim(0,50)
pΜέγεθος κελιού (Cell size (m))
CellArea = res(density1) %>% prod() # Εναλλακτικός τρόπος υπολογισμού του cell size (m)
sqrt(CellArea) # μέγεθος πλευράς κελιού[1] 195.9745
Πόσα κελιά έχουν τιμή > 0
[1] 55
Ποιες είναι οι τιμές των κελιών >0
[1] 1 1 1 4 1 2 2 1 12 1 1 36 65 2 3 4 1 2 2 62 1 1 1 1 6 2 1 9 19 2 8 1 1 4 2 1 2 2 4 1 1 3 5 1 1 3 5 2 2 2 1 2 1 1 1
Τι εμβαδόν καλύπτουν αυτές οι 304 εμπορικές επιχειρήσεις που βρίσκονται σε 55 κελιά?
[1] 2112330
Επαναλαμβάνουμε όλα τα παραπάνω βήματα ξανά, για τη περιοχή της Χίου.
shops_all = osm.getPOI_usingbb(bb2, inkey ="shop" )
shops_chios = osm.combineShops(shops_all)
mypol_bb = osm.osmdata_result_2_bbox_pol(shops_all) %>% st_transform("EPSG:4326") # πολύγωνο Χίου
pol_chios = sf::st_intersection( mypol_bb, acto )
r_chios = rast( vect(pol_chios), res=size ) %>% raster() # Άδειο raster για χρήση στη συνέχεια.Μετατροπή όλων των δεδομένων της Χίου σε wgs:2100
r_chios = projectRaster(r_chios, crs = "EPSG:2100")
shops_chios = shops_chios %>% st_transform("EPSG:2100")
pol_chios = pol_chios %>% st_transform("EPSG:2100")Δεδομένα μέχρι στιγμής:
#par(mar = rep(0, 4)) # Remove all margins
plot(st_geometry(pol_chios))
plot(st_geometry(shops_chios), add=T)Πίνακας συχνοτήτων των επιχειρήσεων στη περιοχή μελέτης.
# A tibble: 32 × 3
shop n freq
* <chr> <int> <dbl>
1 <NA> 20 22.7
2 clothes 11 12.5
3 supermarket 11 12.5
4 bakery 4 4.55
5 mobile_phone 4 4.55
6 computer 3 3.41
7 convenience 3 3.41
8 electronics 3 3.41
9 hairdresser 2 2.27
10 health_food 2 2.27
# ℹ 22 more rows
# A tibble: 32 × 3
shop n freq
* <chr> <int> <dbl>
1 <NA> 20 22.7
2 clothes 11 12.5
3 supermarket 11 12.5
4 bakery 4 4.55
5 mobile_phone 4 4.55
6 computer 3 3.41
7 convenience 3 3.41
8 electronics 3 3.41
9 hairdresser 2 2.27
10 health_food 2 2.27
11 shoes 2 2.27
12 sports 2 2.27
13 travel_agency 2 2.27
14 bag 1 1.14
15 books 1 1.14
16 butcher 1 1.14
17 car 1 1.14
18 car_repair 1 1.14
19 deli 1 1.14
20 florist 1 1.14
21 gift 1 1.14
22 greengrocer 1 1.14
23 jewelry 1 1.14
24 kiosk 1 1.14
25 motorcycle 1 1.14
26 optician 1 1.14
27 paint 1 1.14
28 pastry 1 1.14
29 pet 1 1.14
30 telecommunication 1 1.14
# ℹ 2 more rows
Επιλογή μόνο των εμπορικών επιχειρήσεων
Raster με πλήθος εμπορικών επιχειρήσεων αν κελί
Μέγεθος κελιού (Cell size (m))
CellArea2 = res(density2) %>% prod() # Εναλλακτικός τρόπος υπολογισμού του cell size (m)
sqrt(CellArea2) # μέγεθος πλευράς κελιού[1] 197.104
Οπτικοποίηση
1 2 3 7 14
11 5 1 1 1
hist(density2, xlim=c(0,50), main="Ιστόγραμμα συχνοτήτων (Χίος)\n(πλήθος επιχειρήσεων ανά κελί)", ylab="Συχνότητα", xlab="Επιχειρήσεις")Κάθετη συνένωση των δυο data.frame των δυο περιοχών (Μυτιλήνη, Χίος)
# A tibble: 6 × 2
layer Perioxi
<dbl> <chr>
1 NA Mytilini
2 NA Mytilini
3 NA Mytilini
4 NA Mytilini
5 NA Mytilini
6 NA Mytilini
Στατιστικά κελιών ανα περιοχή
data_all %>% group_by(Perioxi) %>%
summarise(kelia_plithosa=n(),
MO=mean(layer, na.rm=T),
kelia_more_1 = sum( layer>1, na.rm=T ) ,
kelia_more_2 = sum( layer>2, na.rm=T ) ,
kelia_more_3 = sum( layer>3, na.rm=T )
)# A tibble: 2 × 6
Perioxi kelia_plithosa MO kelia_more_1 kelia_more_2 kelia_more_3
<chr> <int> <dbl> <int> <int> <int>
1 Chios 1173 2.37 8 3 2
2 Mytilini 500 5.53 31 17 14
Πόσα κελιά έχει το κάθε raster της κάθε περιοχής
Chios Mytilini
1173 500
Ιστογραμμα 1
p = data_all %>%
ggplot( aes(x=layer, fill=Perioxi)) +
geom_histogram( alpha=0.8, position = 'identity') +
xlim(0,30)
pΙστογραμμα 2
library(ggpubr)
gghistogram(
data_all, x = "layer", alpha=0.6,
add = "mean", rug = TRUE,
fill = "Perioxi", palette = c("blue", "green"), add_density = TRUE
)Συγκρίνετε δύο άλλες πόλεις (πχ. Λάρισα, Χανιά, Πάτρα, Ρόδος), επιλέγοντας τις εμπορικές επιχειρήσεις και απαντώντας τα παρακάτω ερωτήματα ανά πόλη: