1 Στόχοι της άσκησης

2 Εισαγωγή

Εντοπισμός της εμπορικής περιοχής σε μια πόλη χρησιμοποιώντας δεδομένα OpenStreetMap. Θα μελετήσουμε το πλήθος, το είδος επιχειρήσεων και τη χωρική κατανομή τους, σε μια υπο-περιοχή της πόλης της Μυτιλήνης, στο Δήμο Μυτιλήνης, στη Λέσβο.

3 Βιβλιοθήκες

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)

4 Μυτιλήνη

4.1 Δεδομένα Μυτιλήνης

Διαβάζουμε το πολύγωνο της ακτογραμμής της Λέσβου. Στη συνέχεια ορίζουμε ένα δικό μας τετράγωνο περίγραμμα (Bounding Box) το οποίο θα χρησιμοποιήσουμε για να κατεβάσουμε δεδομένα.

acto = sf::read_sf("data/Greece_wgs84.shp")  %>% st_transform("EPSG:4326")

Διαθέσημες πόλεις για ανάλυση (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)

size = 0.002

Χρησιμοποιώντας ένα Bounding-Box (ΒΒ) κατεβάζουμε όλες τις επιχειρήσεις από το OpenStreetMap. (Points, Polygons, Multipolygons) οι οποίες βρίσκονται εντός του Bounding-Box. Στη συνέχεια συνενώνουμε όλα τα είδη των επιχειρήσεων σε ένα σημειακό σύνολο-δεδομένων.

shops_all = osm.getPOI_usingbb(bb1, inkey ="shop" )
shops = osm.combineShops(shops_all)
shops
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 )

Δεδομένα μέχρι στιγμής:

#par(mar = rep(0, 4)) # Remove all margins
plot(st_geometry(pol))
plot(st_geometry(shops), add=T)

Κατασκευάζουμε ένα άδειο raster για χρήση στη συνέχεια.

r = rast( vect(pol), res=size ) %>% raster()
r
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

r = projectRaster(r, crs = "EPSG:2100")
shops = shops %>% st_transform("EPSG:2100")
pol= pol %>% st_transform("EPSG:2100")

4.2 Ανάλυση Μυτιλήνης

Πίνακας συχνοτήτων των επιχειρήσεων στη περιοχή μελέτης.

freq1 = osm.getFrequency(shops, inword = "shop", removeNA = F)
freq1
# 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
print(freq1, n=30) # εκτυπωση 30 γραμμών
# 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,]

Συχνότητα εμπορικών επιχειρήσεων

table(shops_ena$shop) %>% as.data.frame() %>% arrange(desc(Freq))
           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

Συνολικό πλήθος εμπορικών επιχειρήσεων

nrow(shops_ena)
[1] 304

Raster με πλήθος εμπορικών επιχειρήσεων αν κελί

density1 = raster::rasterize(shops_ena, r, field=1, fun=sum) # sum

Raster με Ύπαρξη/Απουσία επιχειρήσεων ανά κελί

#presense1 = raster::rasterize(shops_ena, r, field=1) # presense

Μεθοδος για τον υπολογισμό του πινακα αποστάσεων μεταξύ των σημείων. Ειτε για ολα τα σημεια, ειτε για μια κατηγορία σημειων μονο.

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)
}

Ποιες ειναι οι αποστάσεις μεταξύ ολων των εμπορικών επιχειρησεων γενικά ?

myDist(shops_ena )%>% as.vector( ) %>% summary()
   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 

Οπτικοποιηση

plot(st_geometry(pol))
plot(st_geometry(shops_ena), add=T )
plot(density1, add=T)

Πίνακας συχνοτήτων κελιών

table(density1[])

 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

length(density1[density1>0])
[1] 55

Ποιες είναι οι τιμές των κελιών >0

density1[density1>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 κελιά?

length(density1[density1>0]) * CellArea
[1] 2112330

5 Χίος

5.1 Δεδομένα Χίος

Επαναλαμβάνουμε όλα τα παραπάνω βήματα ξανά, για τη περιοχή της Χίου.

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)

5.2 Ανάλυση Χίος

Πίνακας συχνοτήτων των επιχειρήσεων στη περιοχή μελέτης.

freq2 = osm.getFrequency(shops_chios, inword = "shop", removeNA = F)
freq2
# 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
print(freq2, n=30) # εκτύπωση 30 γραμμών
# 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

Επιλογή μόνο των εμπορικών επιχειρήσεων

shops_duo = shops_chios[shops_chios$shop %in% emporikes_epixeiriseis,]

Raster με πλήθος εμπορικών επιχειρήσεων αν κελί

density2 = raster::rasterize(shops_duo, r_chios, field=1, fun=sum) # sum

Μέγεθος κελιού (Cell size (m))

CellArea2 = res(density2) %>% prod() # Εναλλακτικός τρόπος υπολογισμού του cell size (m)
sqrt(CellArea2) # μέγεθος πλευράς κελιού
[1] 197.104

Οπτικοποίηση

plot(st_geometry(pol_chios))
plot(st_geometry(shops_duo), add=T, cex=0.5 )
plot(density2, add=T)

table(density2[])

 1  2  3  7 14 
11  5  1  1  1 
hist(density2, xlim=c(0,50), main="Ιστόγραμμα συχνοτήτων (Χίος)\n(πλήθος επιχειρήσεων ανά κελί)", ylab="Συχνότητα", xlab="Επιχειρήσεις")

data_hist2 = as.data.frame(density2) # Χίος
data_hist2$Perioxi = "Chios"

6 Σύγκριση περιοχών

Κάθετη συνένωση των δυο data.frame των δυο περιοχών (Μυτιλήνη, Χίος)

data_all = rbind(data_hist, data_hist2) %>% as_tibble()
head(data_all)
# 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 της κάθε περιοχής

table(data_all$Perioxi)

   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
)

7 Ερωτήματα

Συγκρίνετε δύο άλλες πόλεις (πχ. Λάρισα, Χανιά, Πάτρα, Ρόδος), επιλέγοντας τις εμπορικές επιχειρήσεις και απαντώντας τα παρακάτω ερωτήματα ανά πόλη:

  1. Πόσα κελιά έχουν τιμή > 0
  2. Μέγεθος κελιού (μέτρα)
  3. Ποια είναι η πιο συχνή εμπορική επιχείρηση σε κάθε πόλη? Με τι συχνότητα και σχετική συχνότητα?
  4. Συνολικό πλήθος εμπορικών επιχειρήσεων σε κάθε πόλη?
  5. Πλήθος κελιών εμπορικών επιχειρήσεων σε κάθε πόλη?
  6. Συνολικό Εμβαδό κελιών (\(μέτρα^2\)) εμπορικών επιχειρήσεων σε κάθε πόλη?
  7. Ποιά είναι η μέση αποστάση μεταξύ ολων των εμπορικών επιχειρησεων γενικά? Τι συμπεραίνετε για κάθε πόλη?