Bunch of Ideas, a little bit of Maths.

I am mathematician doing market research and data science

14 minute read

Tulisan ini masih menjadi satu rangkaian pengolahan data geomarketing yang saya sedang lakukan di kantor beberapa minggu belakangan ini.

Kali ini saya hendak menuliskan tentang apa saja yang saya oprek dari format data geospasial berekstensi .kml dan .kmz.

Sebagai contoh pada tulisan ini, saya akan menggunakan data Google Mymaps yang saya miliki. Saya memiliki data rekapan rute yang dijalankan oleh beberapa orang interviewer saat melakukan survey di kota Solo. Berikut adalah contoh datanya saat ditampilkan di Google Mymaps.

Datanya sendiri disimpan dalam nama Solo Full.kmz. Jika dilihat dari peta di atas, Solo Full.kmz memiliki beberapa layers rute.

Oke, kita mulai ngopreknya…!

Pertama-tama, kita panggil libraries yang terlibat:

library(dplyr) # data wrangling
library(tidyr) # data wrangling
library(sf)    # geospasial
library(mapview) # membuat maps

Langkah pertama adalah dengan mengecek ada berapa layers data pada file .kmz tersebut.

file     = "Solo Full.kmz"
layers   = st_layers(file)

layers

Driver: LIBKML 
Available layers:
                        layer_name geometry_type features fields crs_name
1  Surakarta-Banjarsari-Banyuanyar                      4     11   WGS 84
2    Surakarta-Banjarsari-Kadipiro                      4     11   WGS 84
3    Surakarta-Banjarsari-Keprabon                      6     11   WGS 84
4     Surakarta-Banjarsari-Manahan                      4     11   WGS 84
5     Surakarta-Banjarsari-Nusukan                      6     11   WGS 84
6     Surakarta-Banjarsari-Timuran                      5     11   WGS 84
7          Surakarta-Laweyan-Jajar                      4     11   WGS 84
8     Surakarta-Laweyan-Karangasem                      4     11   WGS 84
9       Surakarta-Laweyan-Sondakan                      4     11   WGS 84
10     Surakarta-Laweyan-Sriwedari                      4     11   WGS 84

Kita dapatkan ada 10 layers yang merupakan 10 rute yang saya miliki di Solo. Selanjutnya saya hendak melihat detail salah satu layer yang ada. Misalkan layer bernama Surakarta-Banjarsari-Banyuanyar.

df_geo = st_read(file,layer = layers$name[1])

Reading layer `Surakarta-Banjarsari-Banyuanyar' from data source 
  `/root/ikanx101.github.io/_posts/geo marketing/post 7/Solo Full.kmz' 
  using driver `LIBKML'
Simple feature collection with 4 features and 11 fields
Geometry type: GEOMETRY
Dimension:     XYZ
Bounding box:  xmin: 110.8022 ymin: -7.542342 xmax: 110.8096 ymax: -7.536938
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

Berikut ini adalah detail object data yang sudah di-read:

df_geo

Simple feature collection with 4 features and 11 fields
Geometry type: GEOMETRY
Dimension:     XYZ
Bounding box:  xmin: 110.8022 ymin: -7.542342 xmax: 110.8096 ymax: -7.536938
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
                             Name description timestamp begin  end altitudeMode
1 Surakarta-Banjarsari-Banyuanyar        <NA>      <NA>  <NA> <NA>         <NA>
2                    SOTO AMPIRAN        <NA>      <NA>  <NA> <NA>         <NA>
3           Sosis Bakar Jaman Now        <NA>      <NA>  <NA> <NA>         <NA>
4          Rumah Makan Sinar Pagi        <NA>      <NA>  <NA> <NA>         <NA>
  tessellate extrude visibility drawOrder icon                       geometry
1          1       0         -1        NA <NA> LINESTRING Z (110.8022 -7.5...
2         -1       0         -1        NA <NA> POINT Z (110.8022 -7.536938 0)
3         -1       0         -1        NA <NA>  POINT Z (110.808 -7.538976 0)
4         -1       0         -1        NA <NA> POINT Z (110.8096 -7.542342 0)

Terlihat membingungkan ya? Ini versi simpelnya:

df_geo |> glimpse()

Rows: 4
Columns: 12
$ Name         <chr> "Surakarta-Banjarsari-Banyuanyar", "SOTO AMPIRAN", "Sosis…
$ description  <chr> NA, NA, NA, NA
$ timestamp    <dttm> NA, NA, NA, NA
$ begin        <dttm> NA, NA, NA, NA
$ end          <dttm> NA, NA, NA, NA
$ altitudeMode <chr> NA, NA, NA, NA
$ tessellate   <int> 1, -1, -1, -1
$ extrude      <int> 0, 0, 0, 0
$ visibility   <int> -1, -1, -1, -1
$ drawOrder    <int> NA, NA, NA, NA
$ icon         <chr> NA, NA, NA, NA
$ geometry     <GEOMETRY [°]> LINESTRING Z (110.8022 -7.5..., POINT Z (110.8022 -7.5369…

Kita bisa membuat maps di R menggunakan layer tersebut:

mapview(df_geo)

Kita juga bisa memandang data tersebut sebagai struktur data frame.

df_geo |> as.data.frame()

                             Name description timestamp begin  end altitudeMode
1 Surakarta-Banjarsari-Banyuanyar        <NA>      <NA>  <NA> <NA>         <NA>
2                    SOTO AMPIRAN        <NA>      <NA>  <NA> <NA>         <NA>
3           Sosis Bakar Jaman Now        <NA>      <NA>  <NA> <NA>         <NA>
4          Rumah Makan Sinar Pagi        <NA>      <NA>  <NA> <NA>         <NA>
  tessellate extrude visibility drawOrder icon                       geometry
1          1       0         -1        NA <NA> LINESTRING Z (110.8022 -7.5...
2         -1       0         -1        NA <NA> POINT Z (110.8022 -7.536938 0)
3         -1       0         -1        NA <NA>  POINT Z (110.808 -7.538976 0)
4         -1       0         -1        NA <NA> POINT Z (110.8096 -7.542342 0)

Perhatikan bahwa baris pertama adalah rute jalan sedangkan baris kedua sampai akhir merupakan titik tempat.

Kita bisa menghitung berapa jarak dari rute jalan pada baris pertama dengan cara:

st_length(df_geo$geometry[1])

1359.684 [m]

Kita dapatkan jarak sepanjang 1.359,684 meter.

Kita juga bisa mengekstrak titik longlat dari kolom geometri sebagai berikut:

df_geo[2:4,] %>% 
  mutate(long = st_coordinates(.)[,1],
         lat = st_coordinates(.)[,2]) |> 
  select(Name,long,lat) |> 
  as.data.frame()

                    Name     long       lat                       geometry
2           SOTO AMPIRAN 110.8022 -7.536938 POINT Z (110.8022 -7.536938 0)
3  Sosis Bakar Jaman Now 110.8080 -7.538976  POINT Z (110.808 -7.538976 0)
4 Rumah Makan Sinar Pagi 110.8096 -7.542342 POINT Z (110.8096 -7.542342 0)

Jika rekan-rekan masih ingat tulisan saya tentang perhitungan jarak dan membuat rute di osrm, kita juga bisa meng-ekspor rute tersebut menjadi format ekstensi .kml sehingga bisa di-upload dan digunakan pada Google Earth atau Google Mymaps.

Sebagai contoh, saya akan gunakan data sebaran rumah sakit di Provinsi Jawa Barat yang saya dapatkan dari situs arcgis Jawa Barat. Berikut adalah contoh datanya:

  kota nama_rs telepon lat long
204 Bogor RS Sentosa 0251 - 7541900 -6.495095 106.7511
248 Bekasi RS Umum Cikarang Medika 021-29568832 -6.251924 107.1451
159 Bekasi RS Umum Medirossa Cikarang 021-89833216 -6.299245 107.1450
267 Kota Bogor RS Ibu dan Anak Sawojajar (0251) 832 4371 -6.588111 106.7937
34 Kota Tasikmalaya RS Umum Islam Hj. Siti Muniroh 0265-323868 -7.378759 108.2328
378 Sukabumi RS Primaya Sukabumi   -6.918458 106.9578
209 Bogor RS Umum Daerah Cileungsi (021) 89934667 -6.432804 107.0487
166 Kota Bekasi RS Umum Masmitra 2184971766 -6.281745 106.9300
337 Kuningan RS Jantung Hasna Medika Kuningan 6.22329E+11 -6.974656 108.4678
360 Indramayu RS Sentra Medika Langut   -6.406972 108.2506

Misalkan saya hendak membuat rute bersepeda yang melalui 10 Rumah Sakit tersebut. Saya bisa menggunakan function osrmTrip() dari library(osrm) dengan input berupa matriks longlat. Berikut caranya:

input = df |> select(long,lat) |> as.matrix()
input

        long       lat
204 106.7511 -6.495095
248 107.1451 -6.251924
159 107.1450 -6.299245
267 106.7937 -6.588111
34  108.2328 -7.378759
378 106.9578 -6.918458
209 107.0487 -6.432805
166 106.9300 -6.281745
337 108.4678 -6.974656
360 108.2506 -6.406972

trips = osrm::osrmTrip(loc = input,osrm.profile = "bike")
trips

[[1]]
[[1]]$trip
Simple feature collection with 10 features and 4 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 106.751 ymin: -7.378762 xmax: 108.4994 ymax: -6.155105
Geodetic CRS:  WGS 84
   start end  duration distance                       geometry
1    204 166 166.08667  39.4165 LINESTRING (106.751 -6.4949...
2    166 209 115.29167  29.5787 LINESTRING (106.93 -6.28175...
3    209 159  87.08500  22.7214 LINESTRING (107.0487 -6.432...
4    159 248  25.33833   6.3237 LINESTRING (107.145 -6.2990...
5    248 360 616.81833 161.0223 LINESTRING (107.1452 -6.251...
6    360 337 381.56667  98.0548 LINESTRING (108.2495 -6.406...
7    337  34 333.86333  85.9802 LINESTRING (108.468 -6.9747...
8     34 378 875.51833 223.6584 LINESTRING (108.2328 -7.378...
9    378 267 297.03167  72.8262 LINESTRING (106.9578 -6.918...
10   267 204  63.42333  15.1833 LINESTRING (106.7937 -6.588...

[[1]]$summary
[[1]]$summary$duration
[1] 2962.023

[[1]]$summary$distance
[1] 754.7655

Saya dapatkan total durasi perjalanan selama 2962 menit dengan total jarak sebesar 754.8 km.

Berikut adalah petanya:

mytrip = trips[[1]]$trip

mapview(mytrip)

Kita bisa save rute tersebut menjadi file .kml agar bisa dimasukkan ke dalam Google Earth atau Google Mymaps, caranya:

titik = 
  df |> 
  st_as_sf(coords = c("long", "lat"), crs = "+proj=longlat +datum=WGS84")

st_write(obj = titik, dsn = "tes rute.kml", driver = "libkml",layer="pt")
st_write(obj = mytrip, dsn = "tes rute.kml", driver = "libkml")

Berikut adalah screenshot setelah saya upload ke Google Mymaps:

if you find this article helpful, support this blog by clicking the ads.

You May Also Enjoy