A naive approach

Calculations

library(shape)
x =c(99.030269,99.03267,99.03266,99.03268,99.03294)
y =c(18.79611,18.79606,18.79609,18.79608,18.79609)
plot(x,y,main="Location of Palm Tree",
     xlab="Longitude",ylab="Lattitude",pch=19)
points(mean(x),mean(y),pch=17,col="red")
plotcircle(mid=c(mean(x),mean(y)),r=0.001,
           col="#99000033",from=0,to=2*pi)
points(99.032669,18.796070,col="blue",pch=17)

legend(99.0305,18.79608,c("Data points","Average point","Googlemap point"),fill=c("black","red","blue"))

The averages:

Lat 18.796086 \(\pm\) 1.8165902^{-5}

Lng 99.0322438 \(\pm\) 0.0011101

Error in meters equivalent: 5.6655772

Mapping GPS data in R

STEP 1: Load the special libraries

leaflet:

Javascript drawing and mapping functions using various map tile sets

library(leaflet)

rgdal: Bindings for the ‘Geospatial’ Data Abstraction Library. Mathematics for 2D/ 3D maps of earth coordinates

library(rgdal)
## Loading required package: sp
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: C:/Users/VV-106/AppData/Local/R/win-library/4.3/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: C:/Users/VV-106/AppData/Local/R/win-library/4.3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.0-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

dplyr: data table manipulation

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

XML: xml parser

library(XML)

HTML tools and widgets

 library(htmlwidgets)
 library(htmltools)

STEP 2: Establish the bounding box

japan_lon <- 138.129731
japan_lat <- 38.0615855

STEP 3: Load in a GIS dataset and filter the data

This is a public dataset of all earthquakes that have occurred in and around Japan between 2001 and 2018. In this example we only plot the strongest earthquakes (those greater in magnitude that 6.0)

data = read.csv("JapanEarthquakes2001-2018.csv")
strongquakes <- data %>% filter(mag >= 6)

STEP 4: Establish plotting color palette

pcolors = rev(rainbow(10,alpha=0.9,s=1))

เราเขียนภาษาไทย

STEP 5: Set up the Map title layer

tag.map.title <- tags$style(HTML("
  .leaflet-control.map-title { 
    transform: translate(-50%,20%);
    position: fixed !important;
    left: 50%;
    text-align: center;
    padding-left: 10px; 
    padding-right: 10px; 
    background: rgba(255,255,255,0.5);
    color: black;
    font-weight: bold;
    font-size: 14px;}"))

title <- tags$div(
  tag.map.title, HTML("Earthquakes of Japan 2001-2018")
)

STEP 6: Draw the map with the data points, legend and title

leaflet() %>%
  setView(lng = japan_lon, lat = japan_lat, zoom = 4) %>%
  addProviderTiles("Esri.WorldStreetMap") %>%
  addCircles(data = strongquakes,
           radius = 0.75,
           color = pcolors[strongquakes$mag]) %>% 
  addLegend(position = "bottomright",
      colors = pcolors[c(6,8,10)],
      labels = c(6,8,10), opacity = 1,
      title = "Magnitude") %>%
  addControl(title, position = "topleft", className="map-title")
## Assuming "longitude" and "latitude" are longitude and latitude, respectively

STEP 7: Enhance the map and data focus

knitr::kable(data[1:20,])
time latitude longitude depth mag magType nst gap dmin rms net id updated place type horizontalError depthError magError magNst status locationSource magSource
2018-11-27T14:34:20.900Z 48.3780 154.9620 35.00 4.9 mb NA 92 5.044 0.63 us us1000hx1b 2018-11-27T16:06:33.040Z 269km SSW of Severo-Kuril’sk, Russia earthquake 7.6 1.7 0.036 248 reviewed us us
2018-11-26T23:33:50.630Z 36.0733 139.7830 48.82 4.8 mww NA 113 1.359 1.13 us us1000hwvf 2018-11-27T16:44:22.223Z 3km SSW of Sakai, Japan earthquake 6.0 6.1 0.071 19 reviewed us us
2018-11-26T13:04:02.250Z 38.8576 141.8384 50.56 4.5 mb NA 145 1.286 0.84 us us1000hwmi 2018-11-26T23:52:21.074Z 26km SSE of Ofunato, Japan earthquake 8.4 9.5 0.156 12 reviewed us us
2018-11-26T05:20:16.440Z 50.0727 156.1420 66.34 4.6 mb NA 128 3.191 0.62 us us1000hwkm 2018-11-26T08:13:58.040Z 67km S of Severo-Kuril’sk, Russia earthquake 9.7 7.8 0.045 151 reviewed us us
2018-11-25T09:19:05.010Z 33.9500 134.4942 38.19 4.6 mb NA 104 0.558 0.61 us us1000hwig 2018-11-25T23:24:52.615Z 9km SW of Komatsushima, Japan earthquake 3.4 10.1 0.132 17 reviewed us us
2018-11-25T03:16:46.320Z 48.4158 155.0325 35.00 4.6 mb NA 85 4.993 0.77 us us1000hwaf 2018-11-25T03:44:22.040Z 263km SSW of Severo-Kuril’sk, Russia earthquake 6.6 1.8 0.058 90 reviewed us us
2018-11-23T14:30:14.510Z 37.1821 141.1721 46.76 5.2 mb NA 140 0.750 0.80 us us1000hvsm 2018-11-26T02:54:51.216Z 29km ENE of Iwaki, Japan earthquake 7.2 6.2 0.064 81 reviewed us us
2018-11-23T07:19:51.110Z 29.3424 142.3121 10.00 4.7 mb NA 130 2.241 0.76 us us1000hvm0 2018-11-23T08:14:59.040Z 250km N of Chichi-shima, Japan earthquake 9.4 1.9 0.090 37 reviewed us us
2018-11-20T20:16:02.790Z 44.4524 148.0753 101.46 4.7 mww NA 85 3.950 1.08 us us1000hul5 2018-11-21T19:05:21.040Z 88km S of Kuril’sk, Russia earthquake 7.7 2.1 0.103 9 reviewed us us
2018-11-20T19:09:48.760Z 30.4087 130.0687 123.00 5.5 mww NA 30 2.794 0.89 us us1000hujf 2018-11-26T23:45:15.717Z 96km WSW of Nishinoomote, Japan earthquake 5.9 1.8 0.042 54 reviewed us us
2018-11-20T15:38:31.270Z 42.0009 142.7654 73.55 4.5 mb NA 114 0.292 0.61 us us1000hx47 2018-11-27T21:03:51.040Z 49km SE of Shizunai, Japan earthquake 9.4 8.5 0.270 5 reviewed us us
2018-11-18T23:02:08.410Z 24.1937 125.2046 30.25 4.5 mwr NA 76 2.019 1.14 us us1000htrh 2018-11-19T00:57:49.040Z 67km S of Hirara, Japan earthquake 8.3 5.1 0.060 27 reviewed us us
2018-11-18T01:51:22.160Z 30.7822 141.9762 10.00 4.8 mb NA 124 2.964 0.81 us us1000htd9 2018-11-18T02:24:53.040Z Izu Islands, Japan region earthquake 9.0 1.9 0.095 36 reviewed us us
2018-11-17T07:40:05.930Z 25.3931 141.0321 115.25 4.8 mb NA 108 1.986 1.03 us us1000ht4f 2018-11-17T07:58:15.040Z 74km NNW of Iwo Jima, Japan earthquake 6.2 7.3 0.097 33 reviewed us us
2018-11-16T12:09:51.380Z 26.3407 143.8582 36.24 4.8 mb NA 139 1.675 0.93 us us1000hspr 2018-11-16T12:36:09.040Z 182km ESE of Chichi-shima, Japan earthquake 10.3 7.6 0.062 80 reviewed us us
2018-11-16T02:49:12.760Z 39.6014 141.9370 49.23 4.5 mb NA 137 1.361 0.69 us us1000hsi3 2018-11-19T05:23:47.346Z 4km SSW of Miyako, Japan earthquake 7.3 7.5 0.108 25 reviewed us us
2018-11-14T10:07:30.980Z 42.6765 141.9846 35.00 4.8 mb NA 75 1.092 0.61 us us1000hrjc 2018-11-14T11:03:03.312Z 31km E of Tomakomai, Japan earthquake 5.5 1.9 0.040 193 reviewed us us
2018-11-14T04:09:41.120Z 34.7663 139.9805 105.50 5.0 mb NA 97 1.483 0.71 us us1000hrh4 2018-11-15T22:08:18.906Z 26km SSE of Tateyama, Japan earthquake 6.3 6.8 0.107 28 reviewed us us
2018-11-13T21:09:33.230Z 45.2679 150.4810 54.13 4.8 mb NA 129 5.628 0.75 us us1000hrde 2018-11-23T06:08:26.040Z 203km E of Kuril’sk, Russia earthquake 10.0 5.7 0.059 89 reviewed us us
2018-11-12T01:44:52.130Z 42.9554 139.2679 8.66 4.5 mb NA 94 2.546 0.62 us us1000hude 2018-11-21T01:18:05.040Z 101km W of Iwanai, Japan earthquake 7.1 3.0 0.086 40 reviewed us us

Reading GPX files

Reading and parsing GPX

library(XML)

gpx_parsed <- htmlTreeParse(file = "schiehallion.gpx", useInternalNodes = TRUE)

Extracting coordinates and elevation

coords = xpathSApply(doc = gpx_parsed, path = "//trkpt", fun = xmlAttrs)
elevation = xpathSApply(doc = gpx_parsed, path = "//trkpt/ele", fun = xmlValue)

Create the data frame

df <- data.frame(
  lat = as.numeric(coords["lat", ]),
  lon = as.numeric(coords["lon", ]),
  elevation = as.numeric(elevation)
)

head(df)

Naive Plot

plot(x = df$lon, y = df$lat, type = "l",
     col = "black", lwd = 3,
     xlab = "Longitude", ylab = "Latitude")

Draw the map

library(leaflet)

leaflet() %>%
  addTiles() %>%
  addPolylines(data = df, lat = ~lat, lng = ~lon, color = "#000000", opacity = 0.8, weight = 3)

Hiking trail in Scotland

Color chooser Indicating elevation

get_color <- function(elevation) {
  
  if (elevation < 200) {
    return("blue")
  }
  
  if (elevation < 400) {
    return("cyan")
  }
  
  if (elevation < 600) {
    return("green")
  }
  if (elevation < 800) {
    return("yellow")
  }
  if (elevation < 1000) {
    return("orange")
  }
  return("red")
}

New dataset with the new variable for color

df_color <- df %>%
  rowwise() %>%
  mutate(color = get_color(elevation))

df_color$last_color <- dplyr::lag(df_color$color)

Revised Map

map <- leaflet() %>% addTiles()
for (color in levels(as.factor(df_color$color))) {
  map <- addPolylines(map, lat = ~lat, lng = ~lon, data = df_color[df_color$color == color | df_color$last_color == color, ], color = ~color)
}
map

Another example

gpx_parsed2 <- htmlTreeParse(file = "spain.gpx", useInternalNodes = TRUE)

coords2 <- xpathSApply(doc = gpx_parsed2, path = "//trkpt", fun = xmlAttrs)
elevation2 <- xpathSApply(doc = gpx_parsed2, path = "//trkpt/ele", fun = xmlValue)

df2 <- data.frame(
  lat = as.numeric(coords2["lat", ]),
  lon = as.numeric(coords2["lon", ]),
  elevation = as.numeric(elevation2)
)

plot(x = df2$lon, y = df2$lat, type = "l",
     col = "black", lwd = 3,
     xlab = "Longitude", ylab = "Latitude")

leaflet() %>%
  addTiles() %>%
  addPolylines(data = df2, lat = ~lat, lng = ~lon, color = "#000000", opacity = 0.8, weight = 3)

Creating maps with diffrent tile sets

library(leaflet)
japan_lon <- 138.129731
japan_lat <- 38.0615855
leaflet() %>%
  setView(lng = japan_lon, lat = japan_lat, zoom = 4) %>%
  addProviderTiles("Stamen.Toner") %>%
  addCircles(data = strongquakes,
           radius = 0.75,
           color = pcolors[strongquakes$mag]) %>% 
  addLegend(position = "bottomright",
      colors = pcolors[c(6,8,10)],
      labels = c(6,8,10), opacity = 1,
      title = "Magnitude") %>%
  addControl(title, position = "topleft", className="map-title")
## Assuming "longitude" and "latitude" are longitude and latitude, respectively