## Laboratory exercise 9

rm(list=ls())

library(leaflet)
library(sf)
library(DataExplorer)
library(lubridate)

setwd('E://')

pnt_data <- read.csv('rec.csv',header = FALSE, 
                     sep='',skip=15)

#Time instants reformatting
doy <- yday(pnt_data$V1)
res1 <- ymd_hms(paste(pnt_data$V1,pnt_data$V2, sep=' '))
hour <- hour(res1)
min <- minute(res1)
sec <- second(res1)

timer <- doy + hour/24 + min/(24*60) + sec/(24*3600)

latitude <- (pnt_data$V3 + 12.843697222222223)*6378137*pi/180
longitude <- (pnt_data$V4 - 131.13274444444443)*6378137*pi/180
height <- pnt_data$V5 - 125.1

data <- as.data.frame(cbind(longitude, latitude, height))

summary(data)

plot(timer, latitude, type='l', col = 'red', 
     ylim = c(min(latitude, longitude,height),
              max(latitude, longitude,height)),
     xlab='time [day in 2015]',ylab='errors [m]')
lines(timer, longitude, col= 'dark green')
lines(timer, height, col= 'blue')
legend('topright', legend=c('latitude', 'longitude','height'),
       col=c('red','dark green', 'blue'), lty=1:1, cex=0.8,
       box.lty=0)

plot_histogram(data)
plot_correlation(data, maxcat = 5L)
plot_boxplot(data, by = "height")
boxplot(data, names=c('easting [m]',
                      'northing [m]',
                      'height [m]'),
        ylab='# of occurences', main='GPS only')

# Define the coordinates

lat2 <- -12.844
lon2 <- 131.133

# Create a spatial object from the data frame
lon = lon2
lat = lat2

points <- st_as_sf(data.frame(lon= lon, 
                              lat = lat), coords = c("lon", "lat")) #, 
#crs = 4326)

# Add tiles to the map
map <- leaflet(points) %>%
  addTiles() %>%
  # Add markers for the two points with popup labels
  addMarkers(lng = ~lon, lat = ~lat, 
             popup = ~paste("Point", round(lat), round(lon)))
map