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)

## Warning: package 'leaflet' was built under R version 4.2.3

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

library(rgdal)

## Warning: package 'rgdal' was built under R version 4.2.3

## Loading required package: sp

## Warning: package 'sp' was built under R version 4.2.3

## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0

## 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.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/rbatz/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/rbatz/AppData/Local/R/win-library/4.2/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)

## Warning: package 'dplyr' was built under R version 4.2.3

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

## Warning: package 'XML' was built under R version 4.2.3

HTML tools and widgets

 library(htmlwidgets)

## Warning: package 'htmlwidgets' was built under R version 4.2.3

 library(htmltools)

## Warning: package 'htmltools' was built under R version 4.2.3

STEP 2: Establish the bounding box

japan_lat <- 138.129731
japan_lon <- 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 5: Establish plotting color palette

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

STEP 6: 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 7: Draw the map with the data points, legend and title

leaflet() %>%
  setView(lng = japan_lat, lat = japan_lon, zoom = 4) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  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

Changing Tile source

leaflet() %>%
  setView(lng = japan_lat, lat = japan_lon, 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

List of Tile providers

List of providers

STEP 8: Enhance the map and data focus

knitr::kable(data[1:10,])

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

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)

Color chooser Indicating elevation

New dataset with the new variable for color

elevation = as.numeric(elevation)
df$d_elevation <-  abs(c(0,elevation)-c(elevation,0))[1:length(elevation)]

df_color <- df %>%  
  mutate(color = case_when(
    d_elevation < 1 ~ '#000099cc',
    d_elevation < 2 ~ "#009999cc",
    d_elevation < 4 ~ "#009900cc",
    d_elevation < 8 ~ "#999900cc",
    d_elevation < 16 ~ "#993300cc",
    TRUE ~ "#990000cc"))
  
head(df_color)

table(df_color$color)

## 
## #000099cc #009900cc #009999cc #990000cc #993300cc #999900cc 
##         2         1         4        45        14        14

Revised Map

leaflet() %>% addTiles() %>%
  addPolylines(lat = ~lat, lng = ~lon, data = df_color,color = ~color) %>% addCircleMarkers(lat = ~lat, lng = ~lon, data = df_color,color = ~color,radius=3)

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)

elevation2 = as.numeric(elevation2)
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")

d2_elevation = (c(0,elevation2) - c(elevation2,0))[1:length(elevation2)]

df2_color <- df2 %>%  
  mutate(color = case_when(
    d2_elevation < 1 ~ '#000099cc',
    d2_elevation < 2 ~ "#009999cc",
    d2_elevation < 4 ~ "#009900cc",
    d2_elevation < 8 ~ "#999900cc",
    d2_elevation < 16 ~ "#993300cc",
    TRUE ~ "#990000cc"))
  
table(df_color$color)

## 
## #000099cc #009900cc #009999cc #990000cc #993300cc #999900cc 
##         2         1         4        45        14        14

leaflet() %>%
  addTiles() %>%
  addPolylines(data = df2_color, lat = ~lat, lng = ~lon, color = ~color, opacity = 0.8, weight = 3) %>%
addCircleMarkers(lat = ~lat, lng = ~lon, data = df2_color,color = ~color,radius=5)

Low Res Maps

options("sp_evolution_status"=2)
library(sp)
library(maps)

## Warning: package 'maps' was built under R version 4.2.3

library(mapdata)

## Warning: package 'mapdata' was built under R version 4.2.3

library(maptools)

## Warning: package 'maptools' was built under R version 4.2.3

## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: FALSE

map("worldHires", c("Thailand"),col="#ffeecc",boundary=T,fill=T,mar=c(.05,.05,par("mar")[4],0.1),border="#FFCC66",xlim=c(95,120))
colors = rainbow(8)
dat = data.frame(
  lat=c(7.2048417,13.7250472,16.6918797,
        18.7961058,19.9218007,
          17.368688,15.9388149),
  lng=c(99.7961507,100.3027573,98.4753097,
        99.0326617, 99.7961507,102.71802718,
   103.807445 ),
  city=c("SK","Bkk", "MSok",
         "CNX","CRai","NK",
         "RoiEt"),
  colinx = c(1:7))

points(dat$lng,dat$lat,
       col=colors[dat$colinx],pch=19,cex=1.5)
points(dat$lng,dat$lat,cex=1.5)
text(dat$lng,dat$lat,dat$city,pos=4 )
text(dat$lng,dat$lat,dat$colinx,pos=2)
legend(108,20,c("1 SK","2 Bkk", "3 MSok",
         "4 CNX","5 CRai","6 NK",
         "7 RoiEt"),fill=colors[dat$colinx])

Mapping in R