00-base.Rmd

---
title: "00-base"
author: "Pietari Pöykkö"
date: "`r Sys.Date()`"
output: 
  pdf_document: 
    toc: yes
---
R Markdown file

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Elevation system N2000
Coordinate system ETRS-TM35FIN

```{r kaynnistys, echo=F, results='hide'}
# Tuodaan ohjelman käyttäytymistä määrittävät käyttäjän asetukset ympäristöön
source("USERCONF.R")

# Juurikansion määrittäminen muuttujaan
if (!requireNamespace("here", quietly = TRUE)) {
  stop("Package \"here\" must be installed to execute.", call. = FALSE)
}
here::i_am("00-base.Rmd")
root <- here::here()


# Listaa tarvittavat paketit
pkgs <- c("rmarkdown","ggplot2","ggstatsplot","ggrepel","ggridges","lidR",
          "rlas","gstat","sf","terra","dplyr","stringr","memoise","data.table",
          "yaml","rnaturalearth","rnaturalearthdata","rlang","fitdistrplus",
          "DescTools","vctrs","scales","knitr","remotes","DBI","tidyr",
          "future.apply","odbc","httr2","jsonlite","broom","this.path")

# Asenna paketit tarvittaessa
if (curl::has_internet()) {
  utils::install.packages(pkgs[!(pkgs %in% installed.packages(noCache = T))],
                          dependencies = TRUE)
  if (!("rnaturalearthhires" %in% installed.packages())) {
    remotes::install_github("ropensci/rnaturalearthhires")
  }
}

# Ainoat paketit jotka aktivoidaan kokonaan käytön helpottamiseksi
# (muita käytetään eksplisiittisesti, paketti::funktio)
library("ggplot2"); library("ggridges") # More compact plot code
library("future") # To use %seed%
library("DescTools") # To use range operators (e.g. `%[]%`)

for (pkg in pkgs) {
  tryCatch(loadNamespace(pkg),
           error = function(cond){
             warning("Conflict loading namespace ",pkg," from local library:\n",
                     cond,"Package ",pkg," not loaded from local library")})
}
rm(pkg)


# Aseta muuttujat pääkansioiden poluille
# (tallennetaan yhteen muuttujaan jotta ympäristö pysyy siistinä, D = directory)
D <- list(
  data     = file.path(root, "aineisto"),
  koodi    = file.path(root, "R"),
  tulokset = file.path(root, "output"),
  doc      = file.path(root, "doc"),
  # Usein käytetyt alakansiot
  secrets  = file.path(root, "aineisto", "secrets"),
  c        = file.path(root, "R","cache")
)

# Muut tärkeä muuttujat:
# Koordinaattireferenssisysteemin (crs) koodi, joka käytössä Suomessa
# ja käytettävissä aineistoissa, EPSG:3067 = ETRS-TM35FIN 
crs = 3067 

# Aja funktiot kansiosta ympäristöön
files_sources <- list.files(
  file.path(D$koodi,"funktiot"), full.names = T, recursive = T)
invisible(sapply(files_sources, source)) # Aja funktioiden tiedostot
rm(files_sources) # Siivous

# Tarkistetaan annetut kansiot, jotta perusrakenne olisi oikein
# (Myöhemmin asioita tallennetaan myös näiden alakansioihin, joissa tarkistukset pitää tehdä erikseen)
# (Vaimennetaan listan läpikäynnin turhat tekstit)
invisible(capture.output(lapply(D, check_dir, tilde = root, ensure = T)))

# Tarkistetaan myös vaaditut kirjautumistiedot jne.
check_db_info(D$secrets,"db_info.yaml")
check_key(D$secrets, "MML_avain.txt", verbose = T)

rm(pkgs)
```



```{r asemien-lataaminen-pohjavesitietokannasta}
# Ladataan asemataulu tietokannasta, jos sitä ei löydy levyltä oletetusta
# sijainnista, tai jos käyttäjä määrittänyt niin.
dir_asemataulu <- file.path(D$data, "asemataulu.csv")
if (USERCONF$pakota_lataus_asematiedot || !file.exists(dir_asemataulu)) {
  asemataulu <- hae_asematiedot(savepath = dir_asemataulu, tstamp = TRUE)
  cat("Asematiedot haettu suoraan tietokannasta\n")
  # Jos asemat ladataan uusiksi, tulee myös putket päivittää
  USERCONF$pakota_lataus_putkitiedot <- TRUE 
} else {
  # Muussa tapauksessa luetaan levyllä oleva taulu
  asemataulu <- read.csv(
    dir_asemataulu,
    colClasses = c(asema_id = "numeric",
                   tunnus = "character", nimi = "character"),
    na.strings = c("","NA"),
    fileEncoding = "UTF-8-BOM")
  cat("Asemataulu ladattu tiedostosta levyltä:",dir_asemataulu,"\n")
  USERCONF$pakota_lataus_asematiedot <- FALSE
}
rm(dir_asemataulu)

# Draw all monitoring stations, before dropping ones with no pipes
as_kartta <- piirra_asemakartta(asemataulu,
  file.path(D$tulokset,"all_stations_map.png"),
  main = paste(
    "All groundwater-related hydrological\nmonitoring stations, n =",nrow(asemataulu)))

```

```{r paikkojen-lataaminen-pohjavesitietokannasta}
# Ladataan paikkojen ja putkien tiedot vasta asemien laajuuden laskemisen
# jälkeen, koska tällöin pudotetaan pois asemat joilla ei ole yhtään havaintoputkea

### Tehdään putkille sama kun asemille ladattaessa
dir_putket <- file.path(D$data, "kaikki_putket.csv")
if (USERCONF$pakota_lataus_putkitiedot || !file.exists(dir_putket)) {
  # Haetaan putket kaikille asemataulun asemille
  kaikki_putket <- hae_putket(asemataulu[["asema_id"]], savepath = dir_putket)
  cat("Putkien tiedot haettu suoraan tietokannasta")
  # Päivitetään putkien latauksen lippu, koska putket päivitettiin
  USERCONF$pakota_lataus_putkitiedot <- TRUE
} else {
  kaikki_putket <- read.csv(
    dir_putket,
    colClasses = c(asema_tunnus = "character", asema_nimi = "character",
                   tunnus = "character",
                   asema_id = "numeric", putki_id = "numeric"),
    na.strings = c("","NA"),
    fileEncoding = "UTF-8-BOM")
  cat("Putkien tiedot ladattu tiedostosta levyltä:",dir_putket,"\n")
  USERCONF$pakota_lataus_putkitiedot <- FALSE
}

### Sama paikoille (ei rajattu pelkkiin putkiin)
# 1st load location type table (required for setting factor levels if local file is used)
PaikkaTyyppi_path <- file.path(D$data,"PaikkaTyyppi.csv") # Define path
if (USERCONF$pakota_lataus_paikkatiedot || file.exists(PaikkaTyyppi_path)) {
    PaikkaTyyppi <- read.csv(PaikkaTyyppi_path) # Read from file if found or forced
} else {
  # If file not found, query from db, and save to disk for future use.
  PaikkaTyyppi <- DBI::dbGetQuery(
    .connect_db(D$secrets),
    p(readLines(f.path(D$secrets,"SQL","hae_paikat#PaikkaTyyppi.txt")),#Query from file
      collapse = " "))
  data.table::fwrite(PaikkaTyyppi,
                     file = file.path(PaikkaTyyppi_path),
                     na = NA, bom = T)
}
rm(PaikkaTyyppi_path)


dir_paikat <- file.path(D$data, "kaikki_paikat.csv")
if (USERCONF$pakota_lataus_paikkatiedot || !file.exists(dir_paikat)) {
  # Haetaan putket kaikille asemataulun asemille
  kaikki_paikat <- hae_paikat(asema_id_lista = asemataulu[["asema_id"]],
                              savepath = dir_paikat, tstamp = T)
  cat("Havaintopaikkojen tiedot haettu suoraan tietokannasta")
  # Päivitetään paikkojen latauksen lippu, koska paikat päivitettiin
  USERCONF$pakota_lataus_paikkatiedot <- TRUE
} else {
  kaikki_paikat <- read.csv(
    dir_paikat,
    colClasses = c(paikka_id = "numeric", asema_id = "numeric",
                   asema_tunnus = "character", asema_nimi = "character",
                   tunnus = "character"),
    na.strings = c("","NA"),
    fileEncoding = "UTF-8-BOM")
  
  # Convert the loaded structure type col into factor, as when loading from db.
  kaikki_paikat$PaikkaTyyppi_Id <- as.factor(kaikki_paikat$PaikkaTyyppi_Id)
  levels(kaikki_paikat$PaikkaTyyppi_Id) <- split(PaikkaTyyppi$PaikkaTyyppi_Id,PaikkaTyyppi$Selite)
  
  cat("Havaintopaikkojen tiedot ladattu tiedostosta levyltä:",dir_paikat,"\n")
  USERCONF$pakota_lataus_paikkatiedot <- FALSE
}
rm(dir_putket,dir_paikat)
```

```{r asemien-laajuuden-laskeminen}
# Asemien laajuuden laskenta. Pudottaa asemista myös pois sellaiset joilla ei
# ole yhtään havaintoputkia, koska ensisijaisesti kiinnostuttu niistä.
source(file.path(D$koodi, "laske_asemien_laajuus.R"))

# Jos asema mahtuu <100km2 alueelle, voitaisiin KM2 aineisto ladata asemalle yhtenäisenä.
# (eli laajuuden raja (puolikas sivun mitta) < 5000m)
# KM2 kuitenkin ladataan paikkakohtaisesti, joten laajuudelle ei ole käyttöä


# Draw monitoring stations which have monitoring pipes
piirra_asemakartta(asemataulu, file.path(D$tulokset,"stations_w_pipes_map.png"),
                   main = paste(
                     "Hydrological monitoring stations\nwith pipes, n =",
                     nrow(asemataulu)))
# Add remaining stations to old plot with different color
as_kartta + geom_point(data = data.frame(asemataulu$lon, asemataulu$lat),
               mapping = aes(x = asemataulu$lon, y = asemataulu$lat, color="points"),
               shape = "circle",
               colour = "green4",
               size = 3,
               alpha = 0.6) +
  ggtitle(paste("Hydrological monitoring stations"))
  ggsave(file.path(D$tulokset,"stations_combined_map.png"),
         height = 13, width = 8, units = "cm")
rm(as_kartta) # Cleanup
```



```{r ref-mittausten-tarkkuuden-arviointi}
# Load national reference points (if not yet loaded)
source(file.path(D$koodi,"tuo_kiintopisteet.R"))

# Import raw measured point data
ref_raw <- readxl::read_excel(
  file.path(D$data,"ref_mittaukset/ref_mittaukset.xlsx"),
  sheet = "raw_data", skip = 1)
# Calculate corrected value if actual and reported antenna pole lengths differ
ref_raw[["Elev"]] <-
  ref_raw$Elev_raw + (ref_raw$dev_varsi - ref_raw$tod_varsi)


# Calculate differences to some of the national accurate reference points to 
# ensure correct functioning of the device.
source(file.path(D$koodi, "vertaa_mittauksia_kiintopisteisiin.R"))

```

```{r ref-mittausten-kasittely-ja-DTM-arviointi}

# Load and handle reference measurements from monitoring pipes
source(file.path(D$koodi, "tuo_&_kasittele_ref_putkimittaukset.R"))

# Plot a map on which four stations the reference measurement campaign was conducted
save2(points_on_fin(ref_mittaukset$East, ref_mittaukset$North , main=NULL),
      f.path(D$tulokset,"ref_mittaukset","ref_measurement_map"),wid=8, hei=12)


# Calculate the differences between various DTMs and own GPS measurements to
# evaluate which DTM to use as the main source of ground elevation
source(file.path(D$koodi, "vertaa_DTMia_mittauksiin.R"))

```



```{r ELYjen-kenttamittausten-tuonti-ja-liittaminen}
# Tuodaan ELYjen kenttämittausten manuaalisesti koostettu taulukko.
# HUOM! Kaikilta asemilta ei ole tehty tätä kartoituskierrosta

# Jos löytyy taulun versio, joka on jo kertaalleen ajettu tällä koodilla läpi, käytetään sitä.
# (Voidaan myös pakottaa ajamaan paikka_id kohdennus uudestaan)
if (USERCONF$DEV$pakota_kenttamittausten_id_kohdennus) {
  message("ELYjen kenttämittausten tietoihin haetaan paikka ID:t uudelleen käyttäjän määrittelemänä")
  source(file.path(D$koodi, "kirjaa_paikkaidt_kenttamittauksiin.R"))
  # Luo tiedoston seuraavaa vaatimusta varten, eli ei pitäisi koskaan
  # vahingossa ajaa tätä kahdesti (alla oleva else)
}

if (file.exists(file.path(D$data, "kenttamittaukset+idt.csv"))) {
  kenttamittaukset <- data.table::fread(
  file = file.path(D$data, "kenttamittaukset+idt.csv"),
  header = TRUE,
  na.strings = c("NA", ""), blank.lines.skip = TRUE,
  keepLeadingZeros = TRUE,
  encoding = "UTF-8")
} else {
  # Jos valmista tiedostoa ei löydy, otetaan raakaversio josta paikka_id:t puuttuu,
  # ja ajetaan skriptin läpi, joka luo koodia varten muotoillun version.
  # Koodi lisää kenttämittauksiin paikka_id:t yhdistelemällä tietokannan tietoja.
  # Tätä käytetään lähinnä kehitysvaiheessa, jotta saadaan yhdistettyä
  # havaintopaikkojen paikka_id:t ELYjen mittauksiin, ja havaitaan mahd ongelmat.
  # Kun id:t on kaikilla paikoille joilta on kenttämittauksia, ne tallennetaan
  # yllä if-kohdassa luettavaan kenttamittaukset+idt.csv-tiedostoon
  source(file.path(D$koodi, "kirjaa_paikkaidt_kenttamittauksiin.R"))
}


```

```{r pudota-huonot-ELY-mittaukset}
# If certain measurements done by ELYs are not deemed to correct, they are
# discarded here. This is based on background investigations outside of this
# project, and are coded manually.
#Yli 50m muutos jota ei vielä ole viety tietokantaan, ja siirtymän ei ole todettu olevan OK.
huonot <- c(9985, 7836, 7851, 7856,
            7811, 56597, 43975, 7727, 7725, 43463)
# Drop bad
kenttamittaukset <- kenttamittaukset[!kenttamittaukset$paikka_id %in% huonot, ]
rm(huonot)
```

```{r ELY-mittausten-taydennys-omilla-mittauksilla}
#### Supplement ELY measurements with own, known accurate measurements.

# First, check differences between own and ELY measurements
# (Visualize differences to own measurements)
comp_ELY <- merge(ref_mittaukset, kenttamittaukset, by = "paikka_id")
comp_ELY[,"tunnus.y"] <- NULL
comp_ELY <- dplyr::rename(comp_ELY, tunnus = tunnus.x)

comp_ELY <- comp_ELY[!grepl("AM",comp_ELY$tunnus),] #Drop autom.dev. meas. (some pipe lens duped)

# Compare msrmnts of pipe visible section on mossy ground.
hist(comp_ELY[comp_ELY$mittausalusta == "moss", "putki_ely"] - 
       (comp_ELY[comp_ELY$mittausalusta == "moss", "putki_ref_raw"]),
      main = NA, breaks = 20,
     xlab = "ELY pipe vis.len. – ref. meas. pipe vis.len. (mossy surface) [m]")
# The negative differences have been confirmed to result from measuring the height
# from the protective pipe rim instead of the pipe proper.

# Drop duplicated msrmnts of the same locations for the following comparisons.
comp_ELY_loc <- comp_ELY[!(duplicated(comp_ELY$East) & duplicated(comp_ELY$North)), ]
# Drop a few of the most inaccurate ref.measurements
comp_ELY_loc <- comp_ELY_loc[comp_ELY_loc$Vt_Prec < 0.09, ]

# Compare measured ground elevation values
hist(comp_ELY_loc$kor_ely - comp_ELY_loc$Elev, breaks = 30,
     xlab = "ELY elev. – ref. meas. elev. [m]", main = NA)
# Compare measured position coordinates
hist(terra::distance(
  cbind(comp_ELY_loc$East, comp_ELY_loc$North),# own meas.
  cbind(as.numeric(comp_ELY_loc$uusi_lon), as.numeric(comp_ELY_loc$uusi_lat)), # ELY meas.
  lonlat = F, pairwise = T), breaks = 30,
  xlab = "Horizontal difference between ref. meas. and ELY meas. [m]", main = NA)



#### Overwrite appropriate values with own measurements into yhd_mittaukset

# Init variable (ELY measurements as the base)
yhd_mittaukset <- data.frame(kenttamittaukset, source = "ELY")
# Prep also a col for vertical precision from own measurements
yhd_mittaukset[["Vt_Prec"]] <- NA

# Add places from ref measurements which are not in ELY measurements
# Rename cols to enable bind_rows
ow_refmit <- dplyr::rename(
  ref_mittaukset[ , c("paikka_id","tunnus","East","North","Elev","putki_ref","Vt_Prec")],
  uusi_lon = East, uusi_lat = North,
  kor_ely = Elev, putki_ely = putki_ref)
ow_refmit <- data.frame(ow_refmit, source = "Pietari") # Add col to indicate source

# Add rows for places which are missing in ELY measurements but included in own
yhd_mittaukset <- dplyr::bind_rows(
  yhd_mittaukset, 
  ow_refmit[!ow_refmit$paikka_id %in% kenttamittaukset$paikka_id, ])

# Loop through the own reference measurements to replace ELY measurements
for (id in ow_refmit$paikka_id) {
  # For all columns (except "source")
  for (va in names(ow_refmit)[(names(ow_refmit) != "source")]) {
    # Select own measurement value if exists, otherwise select ELY measurement
    yhd_mittaukset[yhd_mittaukset$paikka_id == id, va] <-
      ifelse(!is.na(ow_refmit[ow_refmit$paikka_id == id, va]),
             ow_refmit[ow_refmit$paikka_id == id, va],
             kenttamittaukset[kenttamittaukset$paikka_id == id, va])
    # Also update source information (no empty spots in own measurements, so
    # updating for all columns is OK, albeit redundant computation)
    if (!is.na(ow_refmit[ow_refmit$paikka_id == id, va]))
      yhd_mittaukset[yhd_mittaukset$paikka_id == id, "source"] <- "Pietari"
  }
}
rm(ow_refmit, id, va)

# Rename cols of the combined var to indicate that they are not only ELY meas.
yhd_mittaukset <- dplyr::rename(yhd_mittaukset,
                    yhd_lon = uusi_lon, yhd_lat = uusi_lat,
                    yhd_maa = kor_ely, yhd_putki = putki_ely)

```

```{r yhdista-uudet-mittaukset-kaikkiin-paikkoihin}
# Combine own + ELY measurements with place info. Select database value if missing

# 1s load places from database which are in new measurements but not in orig. fetch
if (USERCONF$pakota_lataus_paikkatiedot == T) {
  lisapaikat <- hae_paikat(
    paikka_id_lista = setdiff(yhd_mittaukset$paikka_id, kaikki_paikat$paikka_id),
    savepath = file.path(D$data,"lisapaikat_yhd_mittauksissa.csv"))
} else {
  lisapaikat <- read.csv(file.path(D$data,"lisapaikat_yhd_mittauksissa.csv"),
                         colClasses = c(asema_tunnus="character",asema_nimi="character",
                              tunnus="character",asema_id="numeric",paikka_id="numeric"),
                         na.strings = c("","NA"),
                         fileEncoding = "UTF-8-BOM")
}

kaikki_paikat <- rbind(kaikki_paikat, lisapaikat)
rm(lisapaikat)

# Rename "source" col prior to joining to indetify what it refers to
yhd_mittaukset_temp <- dplyr::rename(yhd_mittaukset, yhd_source = source)

# Yhdistetään kenttämittausten taulusta relevantit sarakkeet kaikkien paikkojen tauluun
# Jos kenttämittauksissa on mukana paikkoja, joita ei ole kaikki_paikat
# taulussa, ne jäävät yhdistämättä.
# (joillakin asemilla ELYn kartoitus kattaa asemiin linkittämättömiä kohteita)
kaikki_paikat <- dplyr::left_join(kaikki_paikat,
                    yhd_mittaukset_temp[ ,c("yhd_lat","yhd_lon",
                                       "yhd_putki","yhd_maa","Vt_Prec",
                                       "yhd_source","paikka_id")],
                    by = "paikka_id")
rm(yhd_mittaukset_temp)

# Combine the coords into full cols, preferring yhd_ coords, as more precise
kaikki_paikat["comb_lon"] <- ifelse(!is.na(kaikki_paikat$yhd_lon),
                                    kaikki_paikat$yhd_lon, kaikki_paikat$lon)
kaikki_paikat["comb_lat"] <- ifelse(!is.na(kaikki_paikat$yhd_lat),
                                    kaikki_paikat$yhd_lat, kaikki_paikat$lat)


# Calculate a column also for own+ELY measurement pipe tops (easier use later)
kaikki_paikat["yhd_yla"] <- kaikki_paikat$yhd_maa + kaikki_paikat$yhd_putki


```

```{r valmistele-uudet-koordinaatit-tietokantaan-vientiin}
# Perform checks on how much the coordinates would change if uploaded to database
koord_upd_dir <- check_dir(D$tulokset,"koordinaattipaivitykset",ensure=T)
#### EROT KONSULTTIEN PAIKOISSA ############################################
### Check how much coordinates will change on consultant locations
# Get consultant places file
konsulttipath <- file.path(D$data,"konsulttien_havaintopaikat.xlsx")
# Read in excel sheets
konsulttip1 <- readxl::read_excel(konsulttipath, skip = 7, sheet = 1)
konsulttip2 <- readxl::read_excel(konsulttipath, skip = 7, sheet = 2)
konsulttip3 <- readxl::read_excel(konsulttipath, skip = 7, sheet = 3)
konsulttip <- rbind(konsulttip1,konsulttip2) |> rbind(konsulttip3) # Join sheets
# Automatically match paikka_ids
konsulttip[ ,"paikka_id"] <- NA # Prep paikka_id
for (i in 1:nrow(konsulttip)) {
  if (is.na(konsulttip[i,"paikka_id"])) {
    id <- kaikki_paikat[kaikki_paikat$tunnus == konsulttip[[i,"Rakenne"]] &
                    kaikki_paikat$asema_tunnus == konsulttip[[i,"Asema-\r\ntunnus"]],
                    "paikka_id"]
    if (rlang::is_empty(id)) print(konsulttip[i,c("Asema-\r\ntunnus","Rakenne")])
    else konsulttip[i,"paikka_id"] <- id
  }
}
rm(konsulttip1,konsulttip2,konsulttip3, konsulttipath, i, id)
# Check match validity
stopifnot(!anyNA(konsulttip$paikka_id), !any(duplicated(konsulttip$paikka_id)))

# Select only those consultant locations which have been measured by ELY or own meas.
cat("Ero ELY mittauksen ja vanhentuneen(?) palvelukuvauksen välillä konsulttikohteissa")
konsulttip <- konsulttip[konsulttip$paikka_id %in% yhd_mittaukset$paikka_id, ]
dist_df <- merge(yhd_mittaukset, konsulttip, by = "paikka_id") |>
  dplyr::mutate(vanha_lon=NULL, vanha_lat=NULL, liejuun_ely=NULL, maalajit=NULL, `Asema-\r\ntunnus`=NULL, Nimi=NULL)
dist_df[["dist"]] <- terra::distance(
  cbind(as.numeric(dist_df$yhd_lon),as.numeric(dist_df$yhd_lat)),
  cbind(as.numeric(dist_df$`ETRS-TM35FIN itä`),as.numeric(dist_df$`ETRS-TM35FIN pohj.`)),
  pairwise = T, lonlat = F)
print(dist_df[dist_df$dist > 10, ])

cat("Ero ELY mittauksen ja tietokannan välillä konsulttikohteissa.")
konsulttip <- kaikki_paikat[kaikki_paikat$paikka_id %in% konsulttip$paikka_id, ]
dist_df <- konsulttip  |>
  dplyr::mutate(vanha_lon=NULL, vanha_lat=NULL, liejuun_ely=NULL, maalajit=NULL, `Asema-\r\ntunnus`=NULL, Nimi=NULL)
dist_df[["dist"]] <- terra::distance(cbind(as.numeric(dist_df$lon),as.numeric(dist_df$lat)),
                cbind(dist_df$yhd_lon, dist_df$yhd_lat),
                pairwise = T, lonlat = F)
print(dist_df[dist_df$dist > 10, ])
rm(dist_df, konsulttip)




#### EROT TIETOKANNAN JA ELYJEN VÄLILLÄ ####################################
# Compare ELY coord msrmnts to db to spot possible bad or mislabeled ELY msrmnts
# Retrieve older snapshot of monitoring points info, in case the coordinates are
# updated during this research
db_koord <- data.table::fread(
  file = file.path(D$data,"vanhat_paikkatiedot_2024_03_19.csv"),
  na.strings = "NA",
  data.table = F,
  keepLeadingZeros = T)
# Limit to ones which have field measurements done
db_koord <- db_koord[
  db_koord$paikka_id %in% yhd_mittaukset$paikka_id, ]


# Sort both (old db values and field measurements) for calculating differences
db_koord <- db_koord[with(db_koord, order(paikka_id)), ]
kent <- yhd_mittaukset[with(yhd_mittaukset, order(paikka_id)),
                       c("paikka_id","asema_nimi","asema_tunnus","tunnus","yhd_lat","yhd_lon") ]
# Drop excess comparative places (which were not included in the snapshot)
kent <- kent[kent$paikka_id %in% db_koord$paikka_id, ]

stopifnot(all.equal(db_koord$paikka_id, kent$paikka_id))

# Calculate distances to cached database coordinates
kent[["coord_diff_to_db"]] <- 
  terra::distance(
  cbind(kent$yhd_lon, kent$yhd_lat),
  cbind(db_koord$lon , db_koord$lat),
  lonlat = F, pairwise = T) |> round(2)


# Combine old coordinates recorded in 2022 database snapshot with current db
# coordinates of places not yet present in the snapshot.
# Eliminates the effect of ELYs updating coordinates right after their own measurements
db_koord <- data.table::fread(
  file = file.path(D$data,"2022_kenttamitt_paikkatiedot.csv"),
  na.strings = c("NA","NULL",""),
  data.table = F,
  keepLeadingZeros = T)
# Clean up and fix colnames
db_koord <- dplyr::rename(db_koord, lat_2022 = KoordErTmPohj, lon_2022 = KoordErTmIta,
                          paikka_id = Paikka_Id, tunnus = Tunnus)
db_koord <- db_koord[, !names(db_koord) %in% "PaikkaTyyppi_Id"]
# Join yhd_mittaus data copy
db_koord <- dplyr::left_join(db_koord, kent, by = "paikka_id")
# Calculate distance to old 2022 coordinates
db_koord[["coord_diff_to_2022"]] <- 
  terra::distance(
  cbind(db_koord$yhd_lon, db_koord$yhd_lat),
  cbind(db_koord$lon_2022 , db_koord$lat_2022),
  lonlat = F, pairwise = T) |> round(2)

# Merge the vector of old coordinates with fuller data
kent <- dplyr::left_join(kent, 
                         db_koord[,c("paikka_id","lat_2022","lon_2022","coord_diff_to_2022")],
                         by = "paikka_id")

# Ensure data.table for filtering
kent <- data.table::as.data.table(kent)
db_koord <- data.table::as.data.table(db_koord)

#### See stats on how great is the difference bween db and ELY measurements
# Grep to drop groundfrost and snow places (r or s followed by 1/2 digits)
cat("MAE of coordinate differences (database vs. ELYs):")
DescTools::MAE(na.omit(kent[!grepl("^.*[rs]\\d{1,2}$",kent$tunnus) ,coord_diff_to_db]),0)
cat("95% quantile of coordinate differences (database vs. ELYs)")
quantile(na.omit(kent[!grepl("^.*[rs]\\d{1,2}$",kent$tunnus) ,coord_diff_to_db]),0.95)
cat("Median of coordinate differences (database vs. ELYs)")
median(na.omit(kent[!grepl("^.*[rs]\\d{1,2}$",kent$tunnus) ,coord_diff_to_db]))
hist(log10(kent[!grepl("^.*[rs]\\d{1,2}$",kent$tunnus) ,coord_diff_to_db]),breaks=50)

# Select only larger differences
over50m <- kent[(coord_diff_to_db>50|is.na(coord_diff_to_db))&(coord_diff_to_2022>50|is.na(coord_diff_to_2022)),][order(asema_tunnus,-coord_diff_to_db)]
# Drop if both distances are NA
over50m <- over50m[!(is.na(coord_diff_to_db) & is.na(coord_diff_to_2022)),]
print(over50m[ , list(paikka_id,asema_nimi,tunnus,coord_diff_to_db)])
# Save to file (used outside of this project to investigate bad places with ELYs)
data.table::fwrite(over50m,
      file = file.path(koord_upd_dir,"ELY+refmittaus_suuri_ero.csv"), bom = T)

rm(kent, db_koord, over50m)



#### SAVE RESULTS (after checks) ################
# Select coordinate info to be written into db (outside of this project)
db_batch <- yhd_mittaukset[ ,c("paikka_id","tunnus","yhd_lat","yhd_lon","source")]
db_batch <- dplyr::mutate(db_batch, yhd_lat=round(yhd_lat), yhd_lon=round(yhd_lon))
# Add +2m to groundfrost and snowstick coords to improve map legibility
db_batch[grepl("^.*[rs]\\d{1,2}$",db_batch$tunnus), "yhd_lon"] <-
  db_batch[grepl("^.*[rs]\\d{1,2}$",db_batch$tunnus),"yhd_lon"] + 2
db_batch[grepl("^.*[rs]\\d{1,2}$",db_batch$tunnus), "yhd_lat"] <-
  db_batch[grepl("^.*[rs]\\d{1,2}$",db_batch$tunnus),"yhd_lat"] + 2
data.table::fwrite(na.omit(db_batch), # Save file
  file.path(koord_upd_dir,"koordinaattipaivitykset_final.csv"), bom = T)
rm(db_batch, koord_upd_dir)

```



```{r LAS-aineiston-katalogien-muodostaminen, results='hide'}

### Create LAS catalog for areas where stations are located (full Finland too heavy)

# Link stations manually to LAS data via map sheets to enable loading only req areas
source(file.path(D$koodi, "link_asem_karttaleht_las_tuotantoal.R"))
# Create the catalog, if access to LAS data.
if (USERCONF$ajo_las_alustalla == T) {
  source(file.path(D$koodi, "muodosta_las_katalogi.R"))
}

# This should not be needed, as the main catalog is small enough
if (USERCONF$DEV$experimental == T) {
      # Muodostetaan pääkatalogista suppeampi, paikkakohtainen katalogi
      # Tämä siksi, koska rajattavat alueet ovat niin kaukana toisistaan
      
      h <- USERCONF$loc_dtm_margin # Use consistent marginal
      
      # Muutetaan koko aineiston katalogi tallentamaan LAS klipit omaan kansioonsa 
      lidR::opt_output_files(ktlg) <- paste0(
        check_dir(D$data, "5p_putkittainen", ensure = T),
        "/{XCENTER}_{YCENTER}_{ID}")
      # Leikataan pääkatalogista alue jokaiselle paikalle marginaalin mukaan
      paikat_ktlg <- lidR::clip_rectangle(
        ktlg,
        kaikki_paikat[, "lon"] - h, kaikki_paikat[, "lat"] - h,
        kaikki_paikat[, "lon"] + h, kaikki_paikat[, "lat"] + h)
      rm(h)
      # Paikkakohtaisten klippien katalogi käsittelee tiedostoja erillisinä, koska
      # eivät muodosta yhtenäistä aluetta
      lidR::opt_independent_files(paikat_ktlg) <- TRUE
      
      # Palautetaan katalogien tulosten tallennuspolku tyhjäksi, koska ei haluta
      # samaan kansioon muita mahdollisia käsittelyjen tuloksia
      lidR::opt_output_files(ktlg) <- ""
      lidR::opt_output_files(paikat_ktlg) <- ""
      lidR::opt_progress(paikat_ktlg) <- FALSE # Quiet
      
      # Linkitetään tallennetut klipit varalta paikka_id:isiin.
      # Palauttaa muuttujan "indeksi_5p", jossa tiedostojen yhteydet paikkoihin.
      source(file.path(D$koodi, "link_paikkaid_las_klipit.R"))
}

```

```{r LAS-katalogien-tarkistus}
if (isTRUE(USERCONF$ajo_las_alustalla)) {
  # Tarkistetaan katalogin eheys
  lidR::las_check(ktlg)
  # Tarkistetaan miltä katalogi näyttää
  lidR::plot(ktlg, chunk = T)
  
  
  if (USERCONF$DEV$experimental == T) {
    lidR::las_check(paikat_ktlg)
    lidR::plot(paikat_ktlg, chunk = T)
  }
}
```

```{r LAS-aineisto-korkeusmalleiksi-ja-korkeuden-lukeminen, message=FALSE}
# Muodostetaan katalogista DTM

if (isTRUE(USERCONF$ajo_las_alustalla)) {
  # Create DTMs (type selection based on results of vertaa_DTMia_mittauksiin.R)
  fut_create_DTMs(ktlg,
                  kaikki_paikat$lon, kaikki_paikat$lat,
                  dtm_type = "krig50",
                  ids = kaikki_paikat$paikka_id,
                  folder = file.path(D$data,"LAS_DTM"),
                  clean = USERCONF$luo_kaikki_5p_dtm,
                  reso = USERCONF$DEV$dtm_luonti_solukoko,
                  margin = USERCONF$loc_dtm_margin)
}
# Read elevations of the created DTMs
# Use the more precisely defined yhd_lon/lat, and only if not available, use db coords
if (!exists("dtm_elev")) {
dtm_elev <- read_DTM_elev(file.path(D$data,"LAS_DTM","krig50"),
                             kaikki_paikat$paikka_id,
                             kaikki_paikat$comb_lon, kaikki_paikat$comb_lat,
                             c_name = "dtm_elev") |>
            dplyr::rename(.data = _, paikka_id = id) # Rename id col
}
# Adjust read DTM elevations based on the previously determined median shift
# from the reference measurements
dtm_elev["dtm_elev"] <- dtm_elev["dtm_elev"] + krig50shift

# Join elevations by paikka_id
if ("dtm_elev" %in% names(kaikki_paikat))
  kaikki_paikat$dtm_elev <- NULL # Wipe col if exists from prev run
# Join (/re-join)
kaikki_paikat <- merge(kaikki_paikat, dtm_elev, by = "paikka_id")

```



```{r KM2-lataaminen-ja-lukeminen, echo=TRUE, message=FALSE}

# Code below downloads KM2 elevation model data from station and observation
# location coordinates.
# KM2 was evaluated in vertaa_DTMia_mittauksiin.R to be less accurate than LAS
# data, and thus is not the primary source to propose new pipe top / ground elevs.
# 
# KM2 elevations will still be used in locations where LAS is not yet available 

# Define full KM2 data storage folder
km2_dir <- check_dir(D$data, "KM2", ensure = T)
# Download KM2 data (if forced)
if (USERCONF$pakota_lataus_KM2 == T) {
  source(file.path(D$koodi,"lataa_taysi_km2_aineisto.R"))
}

# Elevations from KM2 data
km2_elevs <- read_KM2_elev(km2_dir, kaikki_paikat$paikka_id,
                           kaikki_paikat$comb_lon, kaikki_paikat$comb_lat) |>
  dplyr::rename(paikka_id = id, km2_elev = elev)

# Join elevations by paikka_id
if ("km2_elev" %in% names(kaikki_paikat))
  kaikki_paikat$km2_elev <- NULL # Wipe col if exists from prev run
# Join (/re-join)
kaikki_paikat <- merge(kaikki_paikat, km2_elevs, by="paikka_id")
rm(km2_elevs)

# Elevations from KM2 not shifted like LAS-DTM elevations.
# Spread of errors in KM2 vs. reference measurements too large.

message("KM2 elevation model elevations read for all points")
message("Ilmoitukset 'proj_create_from_database:' ei pitäisi haitata")
```

```{r KM2-aineiston-luokittelun-kirjaaminen}
# Haetaan MML:ltä saatu tieto KM2-mallin laadun jakautumisesta karttalehdittäin
KM2_qual <- data.table::fread(
  file.path(D$data, "km2_laatuluokat","km2_laatuluokat.txt"),
  fill = TRUE)
# Uudelleennimetään tärkeimmät sarakkeet
colnames(KM2_qual)[1] <- "karttalehti"
colnames(KM2_qual)[2] <- "KM2_laatu"


# Hyödynnetään KM2:n laadun linkityksessä asemien karttalehtilinkitystä

linkitykset_5p[["KM2_laatu"]] <- rep(NA, nrow(linkitykset_5p))  # Alustetaan sarake

for (laatu_i in 1:length(KM2_qual$karttalehti)) {
  # Etsitään, minkä aseman rivi täsmää laatutiedon karttalehteen
  i <- grep(paste0("*",KM2_qual[laatu_i,"karttalehti"]), linkitykset_5p$karttalehti)
  
  if (rlang::is_empty(i)) next # Ohitetaan tyhjä match
  
  for (n in i) {   # Siltä varalta että grep löytää monta osumaa
    # Jos ei vielä ole laatuluokkaa, korvataan NA
    if (is.na(linkitykset_5p[n,c("KM2_laatu")])) 
      linkitykset_5p[n,c("KM2_laatu")] <- KM2_qual[laatu_i,c("KM2_laatu")]
    else
      # Jos on jo luokka, lisätään solun "perään"
      linkitykset_5p[n,c("KM2_laatu")] <- paste(linkitykset_5p$KM2_laatu[n],
                                        KM2_qual[laatu_i,c("KM2_laatu")])
  }
}; rm(i, n, laatu_i)

#### Käsitellään listatut luokkatiedot poikkeustapausten varalta
####
# Jaetaan kerätyt luokat erilleen (luokat "I" ja "II")
luokat_eritelty <- sapply(linkitykset_5p$KM2_laatu, strsplit, split = " ")
# Tiivistetään uniikkeihin
luokat <- sapply(luokat_eritelty, unique)
# Jos on kahta laatuluokkaa sekaisin, merkitään luokaksi huonompi, eli "II"
sekaisin <- sapply(luokat, function(x) if (length(x) > 1) "II" )

# Ajetaan kahden laatuluokan omaavien tiedot kaikkiin tietoihin
for (i in 1:length(luokat)) {
  luokat[i] <- ifelse(length(sekaisin[[i]]) == 0, luokat[i], sekaisin[i])
}; rm(i)

# Yhdistetään linkitystauluun ei-käsiteltyjen tietojen päälle
linkitykset_5p$KM2_laatu <- unlist(luokat)

# Tarkistetaan vielä tapaukset, joissa asema on kahdella tuotantoalueella
for (id in linkitykset_5p$asema_id) {
  # Katsotaan asemia joita esiintyy 2 rivillä
  if (nrow(linkitykset_5p[linkitykset_5p$asema_id == id, ]) > 1) {
    asema <- linkitykset_5p[linkitykset_5p$asema_id == id, ]
    # Jos rivien laatutiedot poikkeavat, muutetaan laadut "II"
    if (unique(asema$KM2_laatu) > 1) {
      linkitykset_5p[linkitykset_5p$asema_id == id, "KM2_laatu"] <- "II"
    }
  }
}; rm(id, asema)


#### Jaetaan laatuluokat kaikkien havaintopaikkojen taulukkoon
####
kaikki_paikat[["KM2_laatu"]] <- rep(NA, nrow(kaikki_paikat)) # Alustetaan sarake

# Poistetaan kahdennetut aseman rivit (jos asema kahden tuotantoalueen alueella)
trim_linkit <- linkitykset_5p[!duplicated(linkitykset_5p$asema_id), ]

# Kohdennetaan laatutiedot asema_id:n perusteella
for (id in unique(linkitykset_5p$asema_id)) {
  kaikki_paikat[kaikki_paikat$asema_id %in% id, "KM2_laatu"] <-
    trim_linkit[trim_linkit$asema_id %in% id, "KM2_laatu"]
}; rm(id)

rm(luokat_eritelty, sekaisin, luokat, trim_linkit) #Clean

```



```{r koordinaattien-siirron-vaikutus-korkeuteen}
# Examine how the uncertainty of coords affects the reliability of elevation 
# obtained from a DTM from a desired point

# Calculate the differences in coordinates between ELY meas and ref.meas
coord_err <- data.frame(matrix(NA, nrow(comp_ELY_loc), ncol=3,
                               dimnames=list(NULL, c("lon","lat","dist"))))
coord_err$lon <- comp_ELY_loc$East - comp_ELY_loc$uusi_lon
coord_err$lat <- comp_ELY_loc$North - comp_ELY_loc$uusi_lat
coord_err$dist <- terra::distance(
  cbind(comp_ELY_loc$East, comp_ELY_loc$North),# own meas.
  cbind(as.numeric(comp_ELY_loc$uusi_lon),as.numeric(comp_ELY_loc$uusi_lat)),#ELY meas.
  lonlat = F, pairwise = T)

# Try to visualise error dist, no good results
fitdistrplus::descdist(collapse::na_rm(coord_err$lon))
fitdistrplus::descdist(collapse::na_rm(coord_err$lat))
fitdistrplus::descdist(collapse::na_rm(coord_err$dist), boot = 1e3)

gamma_fit <- fitdistrplus::fitdist(collapse::na_rm(coord_err$dist), "gamma")
plot(gamma_fit)
unif_fit <- fitdistrplus::fitdist(collapse::na_rm(coord_err$dist), "unif")
plot(unif_fit)
rm(unif_fit, gamma_fit)


# Plot the histogram of ELY measurements distance to the corresponding ref.msrmnts
png(f.path(D$tulokset,"ref_mittaukset","ELY_dist_to_ref.png"),
    width=18,height=13,units="cm",res=300)
hist(coord_err$dist,
     breaks=seq(0,to=ceiling(max(coord_err$dist,na.rm=T)),by=0.1), 
     xlab="ELY measurement distance to ref. meas. [m]", main=NA)
# Add vertical lines at various points of the distance error distribution
abline(v = quantile(coord_err$dist,0.95,na.rm=T), lty=2,lwd=3,col="red")#95th percentile
abline(v = mean(coord_err$dist,na.rm=T), lty=2,lwd=3,col="orange2") # Mean
abline(v = median(coord_err$dist,na.rm=T), lty=2,lwd=3,col="green3") # Median
dev.off()

# Set the elev reading error radius as the 95th percentile of the distance error
r <- quantile(coord_err$dist, 0.95, na.rm=T)


# Errors distributed fixed length away, along a circle circumference
# (sampled, random direction) (NOT ULTIMATELY USED)
x <- seq(0,1, length.out = nrow(kaikki_paikat)) |>
  sample(size = nrow(kaikki_paikat), replace = T)
# Using 95th percentile of distance error as radius:
lon_err <- r*cos(2*pi*x) # Calculate lon and lat error shifts separately
lat_err <- r*sin(2*pi*x)
ggplot()+geom_point(aes(x=lat_err, y=lon_err)) # Confirm that works as intended
rm(x, lon_err, lat_err)


# Uniform distribution inside a circle (the used coord shift method)
a <- runif(nrow(kaikki_paikat), 0, 1)
b <- runif(nrow(kaikki_paikat), 0, 1)
# Swap a and b individual values such that a[i] < b[i]
for (i in 1:nrow(kaikki_paikat)) {
  if (a[i] > b[i]) {
    temp <- a[i]
    a[i] <- b[i]
    b[i] <- temp
  }
}; rm(temp, i)

# Calculate the points in a circle
lon_err <- b*r*cos(2*pi*a/b)
lat_err <- b*r*sin(2*pi*a/b)
ggplot()+geom_point(aes(x=lat_err, y=lon_err)) #Viz the uniform distr w/in circle


bwid <- 6e-3 # Define bandwidth for density plotting in this section

# Read dtm elevs where coords shifted within a circle w/ r of 95th percentile of
# ELY msrmnt distance error
shifted_dtm_elev_maxerr <- read_DTM_elev(f.path(D$data,"LAS_DTM","krig50"),
                                  kaikki_paikat$paikka_id,
                                  kaikki_paikat$comb_lon + lon_err,
                                  kaikki_paikat$comb_lat + lat_err,
                                  c_name = "dtm_elev") |>
  dplyr::rename(.data = _, paikka_id = id)
# Perform the elimination of systematic error
shifted_dtm_elev_maxerr$dtm_elev <- shifted_dtm_elev_maxerr$dtm_elev + krig50shift
shifted_km2_elev_maxerr <- read_KM2_elev(km2_dir, kaikki_paikat$paikka_id,
                                          kaikki_paikat$comb_lon + lon_err,
                                          kaikki_paikat$comb_lat + lat_err)
# Define xlims to be common across three of the plots
xlims <-
  c(round(min(shifted_km2_elev_maxerr$elev - kaikki_paikat$km2_elev, na.rm=T),1)-0.1,
    round(max(shifted_km2_elev_maxerr$elev - kaikki_paikat$km2_elev, na.rm=T),1)+0.1)

png(f.path(D$tulokset,"korkeudet","elev_err_from_coordshift_r=maxerr.png"),
    width=17, height=11, units="cm", res=300)
plot(density(shifted_km2_elev_maxerr$elev - kaikki_paikat$km2_elev,
             na.rm = T, bw = bwid),
     lwd = 2, col = "dodgerblue", axes=F, xlim = xlims, ylim=c(0,35.1),
     main="Shift radius: 95th percentile error distance",
     xlab="Change in elevation [m]"); axis(1); axis(2)
text(x=-0.3, y=8, paste("n =",
       length(na.omit(shifted_km2_elev_maxerr$elev - kaikki_paikat$km2_elev))),
     col = "dodgerblue")
points(density(shifted_dtm_elev_maxerr$dtm_elev - kaikki_paikat$dtm_elev,
               na.rm = T, bw = bwid),
       type="l", lwd = 2)
text(x=0.4, y=3, paste("n =",
       length(na.omit(shifted_dtm_elev_maxerr$dtm_elev - kaikki_paikat$dtm_elev))))
abline(v = 0, col="gray", lty=2)
abline(v = Quantile(kaikki_paikat$dtm_elev - shifted_dtm_elev_maxerr$dtm_elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="black")
abline(v = Quantile(kaikki_paikat$km2_elev - shifted_km2_elev_maxerr$elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="dodgerblue")
legend(0.5, 15, legend=c("Lidar-DTM","KM2"), col=c("black","dodgerblue"), lty=1,lwd=2)
dev.off()


# Do the same but using the error median and mean values as the shifting distance
r <- mean(coord_err$dist, na.rm=T)
lon_err <- b*r*cos(2*pi*a/b)
lat_err <- b*r*sin(2*pi*a/b)
shifted_dtm_elev_meanerr <- read_DTM_elev(f.path(D$data,"LAS_DTM","krig50"),
                                  kaikki_paikat$paikka_id,
                                  kaikki_paikat$comb_lon + lon_err,
                                  kaikki_paikat$comb_lat + lat_err,
                                  c_name = "dtm_elev") |>
  dplyr::rename(.data = _, paikka_id = id)
# Perform the elimination of systematic error
shifted_dtm_elev_meanerr$dtm_elev <- shifted_dtm_elev_meanerr$dtm_elev + krig50shift
shifted_km2_elev_meanerr <- read_KM2_elev(km2_dir, kaikki_paikat$paikka_id,
                                          kaikki_paikat$comb_lon + lon_err,
                                          kaikki_paikat$comb_lat + lat_err)
png(f.path(D$tulokset,"korkeudet","elev_err_from_coordshift_r=meanerr.png"),
    width=17, height=11, units="cm", res=300)
plot(density(shifted_km2_elev_meanerr$elev - kaikki_paikat$km2_elev,
             na.rm = T, bw = bwid),
     lwd = 2, col = "dodgerblue", axes=F, xlim = xlims, ylim=c(0,35.1),
     main="Shift radius: Mean of error distance",
     xlab="Change in elevation [m]"); axis(2); axis(1)
text(x=-0.3, y=20, paste("n =",
       length(na.omit(shifted_km2_elev_meanerr$elev - kaikki_paikat$km2_elev))),
     col = "dodgerblue")
points(density(shifted_dtm_elev_meanerr$dtm_elev - kaikki_paikat$dtm_elev,
               na.rm = T, bw = bwid),
     lwd = 2, type="l")
text(x=0.4, y=4, paste("n =",
       length(na.omit(shifted_dtm_elev_meanerr$dtm_elev - kaikki_paikat$dtm_elev))))
abline(v = 0, col="gray", lty=2)
abline(v = Quantile(kaikki_paikat$dtm_elev - shifted_dtm_elev_meanerr$dtm_elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="black")
abline(v = Quantile(kaikki_paikat$km2_elev - shifted_km2_elev_meanerr$elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="dodgerblue")
legend(0.5, 25, legend=c("Lidar-DTM","KM2"), col=c("black","dodgerblue"), lty=1,lwd=2)
dev.off()


# Elevation error if the coord shift radius is the median error in distance
r <- median(coord_err$dist, na.rm=T)
lon_err <- b*r*cos(2*pi*a/b)
lat_err <- b*r*sin(2*pi*a/b)
shifted_dtm_elev_mederr <- read_DTM_elev(f.path(D$data,"LAS_DTM","krig50"),
                                  kaikki_paikat$paikka_id,
                                  kaikki_paikat$comb_lon + lon_err,
                                  kaikki_paikat$comb_lat + lat_err,
                                  c_name = "dtm_elev") |>
  dplyr::rename(.data = _, paikka_id = id)
# Perform the elimination of systematic error
shifted_dtm_elev_mederr$dtm_elev <- shifted_dtm_elev_mederr$dtm_elev + krig50shift
shifted_km2_elev_mederr <- read_KM2_elev(km2_dir, kaikki_paikat$paikka_id,
                                         kaikki_paikat$comb_lon + lon_err,
                                         kaikki_paikat$comb_lat + lat_err)
png(f.path(D$tulokset,"korkeudet","elev_err_from_coordshift_r=medianerr.png"),
    width=17, height=11, units="cm", res=300)
plot(density(shifted_km2_elev_mederr$elev - kaikki_paikat$km2_elev,
               na.rm = T, bw = bwid),
     lwd = 2, axes=F, col = "dodgerblue", xlim = xlims,
     main="Shift radius: Median of error distance",
     xlab="Change in elevation [m]"); axis(1);axis(2)
text(x=-0.35, y=24, paste("n =",
       length(na.omit(shifted_km2_elev_mederr$elev - kaikki_paikat$km2_elev))),
     col = "dodgerblue")
points(density(shifted_dtm_elev_mederr$dtm_elev - kaikki_paikat$dtm_elev,
             na.rm = T, bw = bwid),
       type="l", lwd = 2,)
text(x=0.35, y=9, paste("n =",
       length(na.omit(shifted_dtm_elev_mederr$dtm_elev - kaikki_paikat$dtm_elev))))
abline(v = 0, col="gray", lty=2)
abline(v = Quantile(kaikki_paikat$dtm_elev - shifted_dtm_elev_mederr$dtm_elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="black")
abline(v = Quantile(kaikki_paikat$km2_elev - shifted_km2_elev_mederr$elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="dodgerblue")
legend(0.5, 25, legend=c("Lidar-DTM","KM2"), col=c("black","dodgerblue"), lty=1,lwd=2)
dev.off()


# And one more with the shift radius as double the max error of ELY msrmnts.
r <- 2 * max(coord_err$dist, na.rm=T)
lon_err <- b*r*cos(2*pi*a/b)
lat_err <- b*r*sin(2*pi*a/b)
shifted_dtm_elev_dblmaxerr <- read_DTM_elev(f.path(D$data,"LAS_DTM","krig50"),
                                  kaikki_paikat$paikka_id,
                                  kaikki_paikat$comb_lon + lon_err,
                                  kaikki_paikat$comb_lat + lat_err,
                                  c_name = "dtm_elev") |>
  dplyr::rename(.data = _, paikka_id = id)
# Perform the elimination of systematic error
shifted_dtm_elev_dblmaxerr$dtm_elev <- shifted_dtm_elev_dblmaxerr$dtm_elev + krig50shift
shifted_km2_elev_dblmaxerr <- read_KM2_elev(km2_dir, kaikki_paikat$paikka_id,
                                         kaikki_paikat$comb_lon + lon_err,
                                         kaikki_paikat$comb_lat + lat_err)
png(f.path(D$tulokset,"korkeudet","elev_err_from_coordshift_r=dblmaxerr.png"),
    width=17, height=11, units="cm", res=300)
plot(density(shifted_km2_elev_dblmaxerr$elev - kaikki_paikat$km2_elev,
               na.rm = T, bw = bwid),
     lwd = 2, axes=F, col = "dodgerblue", ylim=c(0,5),
     main="Shift radius: 2 * max error distance",
     xlab="Change in elevation [m]"); axis(1);axis(2)
text(x=-0.7, y=2.5, paste("n =",
       length(na.omit(shifted_km2_elev_dblmaxerr$elev - kaikki_paikat$km2_elev))),
     col = "dodgerblue")
points(density(shifted_dtm_elev_dblmaxerr$dtm_elev - kaikki_paikat$dtm_elev,
             na.rm = T, bw = bwid),
       type="l", lwd = 2,)
text(x=0.7, y=2.5, paste("n =",
       length(na.omit(shifted_dtm_elev_dblmaxerr$dtm_elev - kaikki_paikat$dtm_elev))))
abline(v = 0, col="gray", lty=2)
abline(v = Quantile(kaikki_paikat$dtm_elev - shifted_dtm_elev_dblmaxerr$dtm_elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="black")
abline(v = Quantile(kaikki_paikat$km2_elev - shifted_km2_elev_dblmaxerr$elev,
           probs = c(.025,.975), na.rm = T), lty=3, lwd=1.5, col="dodgerblue")
legend(0.6, 5, legend=c("Lidar-DTM","KM2"), col=c("black","dodgerblue"), lty=1,lwd=2)
dev.off()


rm(a, b, r, bwid, xlims, lon_err, lat_err, km2_dir)
```



```{r ylapaiden-laskeminen-ja-kokoaminen}
# Luodaan yhdistelmäsarake putken maanpäällisen osan pituudelle, käyttämällä 
# tietokannasta laskettua tietoa ja kenttämittauksia. Priorisoidaan
# kenttämittausta, omat ELYn edelle. Käytetään tulosta laskemaan oletetut
# putken yläpäät muiden lähteiden ilmoittamasta maanpinnasta

# Calculate the visible pipe section length from database values
kaikki_paikat[["db_putki"]] <- kaikki_paikat$Yla - kaikki_paikat$KorkeusMaa 
# Prep table for collecting pipe top section lengths, ground and pipe top elevs
putki_yhd_temp <- cbind(kaikki_paikat[,c("paikka_id","db_putki","KM2_laatu")],
                        data.frame(
                          matrix(NA, nrow(kaikki_paikat), 2,
                                 dimnames=list(NULL,c("putki","putki_source")))))

# Loop through all places and note the used source of pipe top section
for (i in 1:nrow(kaikki_paikat)) {
  # If ELY/own measurement for pipe is found, use that
  if (!is.na(kaikki_paikat[i,"yhd_putki"])) {
    putki_yhd_temp[i,"putki"] <- kaikki_paikat[i, "yhd_putki"]
    putki_yhd_temp[i,"putki_source"] <- kaikki_paikat[i, "yhd_source"]
  }
  else if (!is.na(kaikki_paikat[i,"db_putki"])) { #If no ELY/own msrmnt, get from db
    # Also check that length is OK(10cm-2m) and not exactly 1m (oft-used placeholder)
    if (kaikki_paikat[i,"db_putki"] > 0.1 &&
        kaikki_paikat[i,"db_putki"] < 2 &&
        kaikki_paikat[i,"db_putki"] != 1.000) {
      putki_yhd_temp[i,"putki"] <- kaikki_paikat[i, "db_putki"]
      putki_yhd_temp[i,"putki_source"] <- "DB" # Note database as the source
    } else { # If length in database not positive or exact 1m, set to 1m and note as missing
      putki_yhd_temp[i,"putki"] <- 1
      putki_yhd_temp[i,"putki_source"] <- "MISSING"
    }
  }
  else { #If no ELY msrmnt and no required db info, set to 1 and note source missing
    putki_yhd_temp[i,"putki"] <- 1
    putki_yhd_temp[i,"putki_source"] <- "MISSING"
  }
}; rm(i)


# Prep pipe ground and top elevation cols to be filled (also include db value cols)
putki_yhd_temp <- 
  cbind(putki_yhd_temp,
        kaikki_paikat["Yla"],
        data.frame(matrix(NA, nrow(kaikki_paikat), 3,
                          dimnames = list(NULL,c("yla","yla_source","yla_verif")))),
        kaikki_paikat["KorkeusMaa"],
        data.frame(matrix(NA, nrow(kaikki_paikat), 3,
                          dimnames = list(NULL,c("maa","maa_source","maa_verif")))))
putki_yhd_temp <- dplyr::rename(putki_yhd_temp, yla_db = Yla)#Rename db pipe top col


# Present the confidence intervals (in meters) of various known errors (used as a
# reference on how much the new value has to differ from database to be justified)

# l_own <- own msrmnt tolerance(+/-) Device reported vt.prec checked for each msrmnt
l_ely <- 0.50 #ELY measurements(+/-) (as per a histogram of a few stations around line 318)
l_dtm <- sqrt(
  (DescTools::Quantile( # DTM error low and high quantiles
  ok_refpoints$krig50_error,
  weights = 1/(ok_refpoints$Vt_Prec^2),#Weighted by vertical precision of ref elev
  probs = c(.025, .975)
  ) - # Adjusted by the shifted amount (krig50 err. weighted median)
  krig50shift) ^ 2 + # Squared to combine the uncertainties correctly.
  # Add uncertainty of the effect of imprecise coordinates when DTM elevs are read
  Quantile(kaikki_paikat$dtm_elev - shifted_dtm_elev_maxerr$dtm_elev,
           probs = c(.025,.975), na.rm = T) ^ 2) # also squared
l_dtm[[1]] <- l_dtm[[1]] * -1 # Lower boundary returned back to negative
l_km2 <- sqrt(
  DescTools::Quantile( # KM2 error, low and high quantiles
  ok_refpoints$km2_error,
  weights = 1/(ok_refpoints$Vt_Prec^2),#Weighted by vertical precision or ref elev
  probs = c(.025, .975)) ^ 2 + #Not shifted,as KM2 elevs were not corrected.
  # Add uncertainty of the effect of imprecise coordinates when KM2 elevs are read
  Quantile(kaikki_paikat$km2_elev - shifted_km2_elev_maxerr$elev,
           probs = c(.025,.975), na.rm = T) ^ 2)
l_km2[[1]] <- l_km2[[1]] * -1 # Lower boundary returned back to negative

# OVERWRITE WITH CONSTANT VALUES DETERMINED EARLIER (Pöykkö, 2024)
# (some coords have improved/changed in db, affecting the calculated coord-caused uncert)
l_dtm <- c("2.5%" = -0.407, "97.5%" = 0.407)
l_km2 <- c("2.5%" = -0.390, "97.5%" = 0.556)

# As the model errors were calculated by "true - model", the same order is
# required when checking whether the difference is within the confidence interval


# Loop through all places for ground elevation and pipe top value.
# Handle pipe tops 1st, as they have the potential to be more correct
# Ground elevs are inferred from pipe tops, such that visible pipe sections are
# always updated to be correct (by shifting the ground) 
for (i in 1:nrow(kaikki_paikat)) {
  # Extract existing vars for easier readability
  yla_db     <- kaikki_paikat[i,"Yla"]
  KorkeusMaa <- kaikki_paikat[i,"KorkeusMaa"]
  yhd_yla    <- kaikki_paikat[i,"yhd_yla"]
  yhd_maa    <- kaikki_paikat[i,"yhd_maa"]
  yhd_source <- kaikki_paikat[i,"yhd_source"]
  Vt_Prec    <- kaikki_paikat[i,"Vt_Prec"]
  dtm_elev   <- kaikki_paikat[i,"dtm_elev"]
  km2_elev   <- kaikki_paikat[i,"km2_elev"]
  putki      <- putki_yhd_temp[i,"putki"]
  putki_source <- putki_yhd_temp[i,"putki_source"]
  # Init new vars so they exist (though should get new value by the end of a loop iter)
  yla <- yla_source <- yla_verif <- NA # New value of the pipe top elevation
  maa <- maa_source <- maa_verif <- NA # New value of the ground elevation
  
  # Check if own measurements are available
  if (yhd_source %in% "Pietari") {
    # Check own pipe top against database (vt.precision of the msrmnt as tolerance)
    if (in_prox(yla_db, yhd_yla, Vt_Prec, na_to_f = T)) {
      # If db already within tolerance of own msrmnt, no new values (write db value)
      yla <- yla_db
      yla_source <- "DB" # Write source (to note as reliable)
      yla_verif <- "Pietari" # Note also the verification source
    } else { # Use own values, if difference beyond limit
      yla <- yhd_yla
      yla_source <- "Pietari"
      yla_verif <- "BAD_DB?"
    }
    # Check own ground elevation against database
    if (in_prox(KorkeusMaa, yhd_maa, Vt_Prec, na_to_f = T)) {
      maa <- KorkeusMaa
      maa_source <- "DB"
      maa_verif <- "Pietari"
    } else {
      maa <- yla - putki
      maa_source <- p0(yla_source,"-",putki_source)
      maa_verif <- "BAD_DB?"
    }
  }
  # Check if ELY measurements are available
  else if (yhd_source %in% "ELY") {
    # If ELY's pipe top within tolerance of database value, do additional checks
    if (in_prox(yla_db, yhd_yla, l_ely, na_to_f = T)) {
      # Check other sources if ELY too close to db (some ELY msrmnts already in db)
      if (in_prox(yla_db, yhd_yla, 0.01, na_to_f = T)) {
        # Check LAS-DTM
        if (!is.na(dtm_elev) &&
            isTRUE((yhd_yla - (dtm_elev+putki)) %[]% l_dtm)) {
          yla <- yla_db
          yla_source <- "DB"
          yla_verif <- p0("ELY&(DTM+",putki_source,")")
        # Check KM2
        } else if (!is.na(km2_elev) &&
                   isTRUE((yhd_yla - (km2_elev+putki)) %[]% l_km2)) {
          yla <- yla_db
          yla_source <- "DB"
          yla_verif <- p0("ELY&(KM2+",putki_source,")")
        } else {
          yla <- yla_db
          yla_source <- "DB"
          yla_verif <- "ELY_IN_DB?"
        }
      #If ELY close to db but ELY results havent been uploaded (not too close)
      } else if (!is.na(dtm_elev) &&
            isTRUE((yhd_yla - (dtm_elev+putki)) %[]% l_dtm)) {
        yla <- yla_db
        yla_source <- "DB"
        yla_verif <- p0("ELY&(DTM+",putki_source,")")
      } else if (!is.na(km2_elev) &&
                   isTRUE((yhd_yla - (km2_elev+putki)) %[]% l_km2)) {
        yla <- yla_db
        yla_source <- "DB"
        yla_verif <- p0("ELY&(KM2+",putki_source,")")
      } else { 
        yla <- yla_db
        yla_source <- "DB"
        yla_verif <- p0("ELY±",l_ely,"m")
      }
    # If too large difference to DB, check top elev also against other available elevs
    } else if (!is.na(yhd_yla)) {
        # 1st, is DTM elev available, and if DTM+pipe w/in dtm-tolerance of ELY's
        if (!is.na(dtm_elev) &&
            isTRUE((yhd_yla - (dtm_elev + putki)) %[]% l_dtm)) {
          yla <- yhd_yla
          yla_source <- "ELY"
          yla_verif <- p0("DTM+",putki_source)
        # If DTM not available, check if within range of KM2+pipe elevation
        } else if (!is.na(km2_elev) &&
                   isTRUE((yhd_yla - (km2_elev + putki)) %[]% l_km2)) {
          yla <- yhd_yla
          yla_source <- "ELY"
          yla_verif <- p0("KM2+",putki_source)
        # If neither source available, note no verification source
        } else if (is.na(dtm_elev) && is.na(km2_elev)) {
          yla <- yhd_yla
          yla_source <- "ELY"
          yla_verif <- "BAD_DB?(NO_VERIF)"
        # "Checks failed" if all other available elevs differ too much from ELY's
        } else {
          yla <- yhd_yla
          yla_source <- "ELY"
          yla_verif <- "CHECKS_FAILED"
        }
    }
  }
  # Check if elev from DTM is available, and that no source was found in prev steps 
  if (!is.na(dtm_elev) && is.na(yla_source)) {
    # 1st check if DTM+pipe is within tolerance of db top value
    if (isTRUE((yla_db - (dtm_elev+putki)) %[]% l_dtm)) {
      yla <- yla_db
      yla_source <- "DB"
      yla_verif <- p0("DTM+",putki_source)
      # If too far, check also against KM2 elevation
    } else if (!is.na(km2_elev) &&
               isTRUE((dtm_elev - km2_elev) %[]% l_km2)) {
      yla <- dtm_elev + putki
      yla_source <- p0("DTM+",putki_source)
      yla_verif <- p0("KM2+",putki_source)
    # If KM2 outside of tolerance, select the DTM elev and note differences
    } else {
      yla <- dtm_elev + putki
      yla_source <- p0("DTM+",putki_source)
      yla_verif <-
        paste("BAD_DB?(KM2_FAIL) DB-pipe∆",round(yla_db - putki - dtm_elev,3),
              ", KM2∆",round(km2_elev - dtm_elev,3))
    }
  }
  # Last option is to use KM2 elevation, if no ground elev found in earlier steps
  if (!is.na(km2_elev) && is.na(yla_source)) {
    if (isTRUE((yla_db - (km2_elev+putki)) %[]% l_km2)) {
      yla <- yla_db
      yla_source <- "DB"
      yla_verif <- p0("KM2+",putki_source)
    } else { # There should not be any more elevation sources left to compare to
      yla <- km2_elev + putki
      yla_source <- p0("KM2+",putki_source)
      yla_verif <- # Check to db failed, but write elev change
        paste("BAD_DB? DB-pipe∆",round(yla_db - putki - km2_elev, 3))
    }
  }
  
  # If none alternative top elev sources are available, use the db value
  # (Should only apply to those which lack the DTM and KM2 files)
  if (all(is.na(c(yhd_yla,yhd_maa,yhd_source, km2_elev,dtm_elev)))) {
    yla <- yla_db
    yla_source <- "DB"
    yla_verif <- "NO_SOURCES"
  }
  
  
  # Get / calculate ground surface elevation
  # Set ground elevations based on pipe top - pipe visible section
  if (is.na(maa_source)) { # Check that no previous ground source (own msrmnts) 
    maa <- yla - putki
    # Set source as top source minus pipe top section
    if (identical(yla_source, putki_source)) {
      maa_source <- yla_source
    } else maa_source <- p0(yla_source,"-",putki_source)
    # Try to additionally verify the ground elev with LAS-DTM and KM2
    if (!is.na(dtm_elev) &&
              isTRUE((maa - dtm_elev) %[]% l_dtm)) {
      maa_verif <- "DTM"
    } else if (!is.na(km2_elev) &&
                     isTRUE((maa - km2_elev) %[]% l_km2)) {
        if (maa_source == "KM2+MISSING-MISSING") {
          maa_verif <- "NO_SOURCES"
        } else {
          maa_verif <- "KM2"
        }
    } else if (is.na(dtm_elev) && is.na(km2_elev)) {
      maa_verif <- "NO_SOURCES"
    } else maa_verif <- "CHECKS_FAILED"
  }
  # Fix source showing as added and then subtracted from the ground source text
  maa_source <- sub("^([A-Z].*)\\+(.*)-\\2$", "\\1", maa_source)
    
  # Assign loop variables to the correct location in the collection table
  putki_yhd_temp[i, "yla"]        <- yla
  putki_yhd_temp[i, "yla_source"] <- yla_source
  putki_yhd_temp[i, "yla_verif"]  <- yla_verif
  putki_yhd_temp[i, "maa"]        <- maa 
  putki_yhd_temp[i, "maa_source"] <- maa_source
  putki_yhd_temp[i, "maa_verif"]  <- maa_verif
}
rm(yla_db,KorkeusMaa, yhd_yla,yhd_maa,yhd_source, Vt_Prec, dtm_elev,km2_elev,
   putki,putki_source, maa,maa_source,maa_verif, yla,yla_source,yla_verif)



# Perform sanity checks on the data
elevtbl_sanitycheck <- function(tbl, skip_yla_tyyppi = NULL) {
  tlb_skip <- tbl[!tbl$PaikkaTyyppi_Id %in% skip_yla_tyyppi, ]
  stopifnot(!anyNA(tlb_skip$yla), !anyNA(tbl$maa)) # Elevs filled
  stopifnot(!anyNA(tlb_skip$yla_source), !anyNA(tbl$maa_source)) #Sources filled
  stopifnot(all(tbl$maa > -5), all(tbl$maa < 1300)) #Ground elevs in [-5;1300]m
  stopifnot(all(tlb_skip$yla > -5),all(tlb_skip$yla < 1300)) #Top elevs [-5;1300]m
  stopifnot(all(tlb_skip$yla - tlb_skip$maa > 0.1)) #Visible sections >10cm
  stopifnot(all(tlb_skip$yla - tlb_skip$maa < 2)) #Visible sections <2m
  stopifnot(all(tlb_skip$yla > tlb_skip$maa)) #All tops above ground
  message("Elevation table sanity check passed")
}
elevtbl_sanitycheck(putki_yhd_temp)#No types skipped, above loop should fill all


# Add identifying information to the table (for possible manual checks and sorting)
putki_yhd_temp <- merge(
  kaikki_paikat[
    ,c("paikka_id","asema_tunnus","asema_nimi","tunnus","PaikkaTyyppi_Id")],
  putki_yhd_temp, by = "paikka_id")

```

``` {r clean-elevation-table}

# Remove top elevs for locations other than wells and monitoring pipes
for (i in 1:nrow(putki_yhd_temp)) {
  # Check that type is NOT monitoring pipe (incl. some moistre meters) or well
  if (!putki_yhd_temp[i, "PaikkaTyyppi_Id"] %in%
      c("Havaintoputki","Kaivo","Maankosteusmittari")) {
    # Set NA to all cols relating to "pipe" top value
    putki_yhd_temp[i,
                   c("db_putki","putki","putki_source",
                     "yla_db","yla","yla_source","yla_verif")] <- NA
  }
}; rm(i) # Clean up

# Drop non-applicable types (no ground elev (and no top elev))
putki_yhd_temp <-
  putki_yhd_temp[!putki_yhd_temp$PaikkaTyyppi_Id %in%
                   c("Pohjavesilammikko","Lähde","Ottamon hana",
                     "Talousvesihana", "Virtaamahavaintopaikka","Luminäyte"), ]
# Drop ground moisture meters (only 25 in the whole database)
putki_yhd_temp <- putki_yhd_temp[!putki_yhd_temp$PaikkaTyyppi_Id %in%
                   c("Maankosteusmittari"), ]



#### Ensure that pairs of manual and autom. pipes have the same elev values
link_am_man <- hae_PaikkaAutomaattiManuaali() # Fetch link table

putki_yhd_temp_backup <- putki_yhd_temp # Back up the table before fixes
# putki_yhd_temp_backup -> putki_yhd_temp

# Pre-emptively set some known exceptions (Hardcoded fixes)
hard_fixes <- tibble::tribble(
  ~val_from, ~val_to,
  60853, 47909, # 0302p5bAM->p5b
  39419, 75430, # (0305)HP3->0305HP3AM
  39422, 75431, # (0305)HP6->0305HP6AM
  39424, 75432, # (0305)HP11->0305HP11AM
  39425, 75440, # (0305)HP12->0305HP12AM
  39423, 75441, # (0305)HP9->0305HP9AM
  39758, 75442, # (0305)HP10->0305HP10AM
  39413, 82160, # (0305)HP1->0305HP1AM
  39414, 82161, # (0305)HP2->0305HP2AM
  39420, 82162, # (0305)HP4->0305HP4AM
  39421, 82163, # (0305)HP5->0305HP5AM
  42041, 85987, # 0801p1b->0801p1bAM
  86023, 9944   # 1209limAM->1209p10(limni)
) |> as.data.frame()
# Apply the hardcode fixes table
for (i in 1:nrow(hard_fixes)) {
  putki_yhd_temp[putki_yhd_temp$paikka_id == hard_fixes[i,"val_to"],
                 c("yla","yla_source","yla_verif","maa","maa_source","maa_verif")] <-
    putki_yhd_temp[putki_yhd_temp$paikka_id == hard_fixes[i,"val_from"],
                   c("yla","yla_source","yla_verif","maa","maa_source","maa_verif")]
  # Unused method (the flexibility is not needed in this case)
  # test <- replace_w_key(test,"paikka_id", toK=hard_fixes[i,"val_to"],
  #   vars=c("yla","maa"),fromK=hard_fixes[i, "val_from"])
}
rm(hard_fixes)


# Set the proritisation order of the available sources (used in the next loop)
yla_src_prio <- c("Pietari","DB","ELY","DTM+ELY","DTM+DB","DTM+MISSING",
                  "KM2+ELY","KM2+DB","KM2+MISSING", NA)
stopifnot(length(unique(putki_yhd_temp$yla_source)) == length(yla_src_prio))
maa_src_prio <- c("Pietari-Pietari","DB-Pietari","DB-ELY","DB","DB-MISSING",
                  "ELY","DTM","KM2")
stopifnot(length(unique(putki_yhd_temp$maa_source)) == length(maa_src_prio))
# Set limit of 1m as the max change between man and AM
diff_lim <- 1 # (stops if exceeded, manual investigations possibly needed)
# Go through the man-AM pairs, select the better found source for both of them
for (i in 1:nrow(link_am_man)) {
  # Get ids from link table
  a_id <- link_am_man[i, "AutomaattiPaikka_Id"]
  m_id <- link_am_man[i, "ManuaaliPaikka_Id"]
  # Get place metadata values (also of those potentially dropped earlier)
  p_a <- putki_yhd_temp[putki_yhd_temp$paikka_id == a_id, ]
  p_m <- putki_yhd_temp[putki_yhd_temp$paikka_id == m_id, ]
  
  suppressWarnings(rm(new_maa, new_yla)) # Ensure fresh vars
  
  # If both am & man have been dropped previously, equality is not checked
  if (nrow(p_a) == 0 && nrow(p_m) == 0) {
    next
  # If top elevs already equal, ensure that ground elevs are equal too
  } else if (isTRUE(all.equal(p_a$yla, p_m$yla))) {
    if (isTRUE(all.equal(p_a$maa, p_m$maa))) {
      next # If both top and ground elevs are equal among man and AM, next 
    } else {
      # Compare ground elevation sources
      for (src in maa_src_prio) {
        # Start comparison once the first prioritised source is matched
        if ((p_a$maa_source == src) || (p_m$maa_source == src)) {
          new_maa <- eq_tester(p_a, p_m, "maa_source",should_eq=src, out_var="maa")
          if (abs(p_a$maa-new_maa) > diff_lim || abs(p_m$maa-new_maa) > diff_lim) stop()
          putki_yhd_temp[putki_yhd_temp$paikka_id == a_id, "maa"] <- new_maa
          putki_yhd_temp[putki_yhd_temp$paikka_id == m_id, "maa"] <- new_maa
          break # Stop looking for further prioritised sources once 1 is found
        }
      }
    }
  } else {
    # If top elevs are not equal, perform checks in the order of prioritised sources
    for (src in yla_src_prio) {
        # Start comparison once the first prioritised source is matched
        if ((p_a$yla_source == src) || (p_m$yla_source == src)) {
          new_yla <- eq_tester(p_a, p_m, "yla_source",should_eq=src, out_var="yla")
          if (abs(p_a$yla-new_yla) > diff_lim || abs(p_m$yla-new_yla) > diff_lim) stop()
          putki_yhd_temp[putki_yhd_temp$paikka_id == a_id, "yla"] <- new_yla
          putki_yhd_temp[putki_yhd_temp$paikka_id == m_id, "yla"] <- new_yla
          break # Stop looking for further prioritised sources once 1 is found
        }
    }
    # Check also ground elevs
    if (isTRUE(all.equal(p_a$maa, p_m$maa))) {
      next
    } else {
      for (src in maa_src_prio) { # Compare ground elev sources
        if ((p_a$maa_source == src) || (p_m$maa_source == src)) {#Compare on 1st matched source
          new_maa <- eq_tester(p_a, p_m, "maa_source",should_eq=src, out_var="maa")
          if (abs(p_a$maa-new_maa) > diff_lim || abs(p_m$maa-new_maa) > diff_lim) stop()
          putki_yhd_temp[putki_yhd_temp$paikka_id == a_id, "maa"] <- new_maa
          putki_yhd_temp[putki_yhd_temp$paikka_id == m_id, "maa"] <- new_maa
          break # Stop looking for further prioritised sources once 1 is found
        }
      }
    }
  }
  
}
rm(a_id,m_id,p_a,p_m, yla_src_prio,maa_src_prio, src, diff_lim)
suppressWarnings(rm(new_maa, new_yla))#might be missing(fresh vars ensured above)


### Set some elevations in the results manually based on more detailed inspections
source(file.path(root,"R","outkorkeus_investigated_fixes.R"))


### Drop some places from the results.
# Create a table copy, where all drops happen
putki_yhd_trim <- putki_yhd_temp
# Create a table where all dropped will be stored
putki_yhd_dropped  <- putki_yhd_temp[0,]

# Do id/station drops. Based on manual investigations. Only affects putki_yhd_trim
# (Some of the now-dropped will be included later in an appendix upload)
source(file.path(root,"R","outkorkeus_investigated_omits.R"))
# Drop snowsticks from the db upload table (sticks missing grnd.elev definition)
putki_yhd_trim <- putki_yhd_trim[!putki_yhd_trim$PaikkaTyyppi_Id %in% "Lumikeppi",]

# All dropped cases above must be handled at a later date by communications and
# manual data updates in cases where updating is deemed to be OK. 


# Code to check that man am pair coords are equal
for (i in seq_len(nrow(link_am_man))) {
  id1 <- link_am_man[i, 1]
  id2 <- link_am_man[i, 2]
  if (!(id1 %in% kaikki_paikat$paikka_id)) {
    warning(paste("Paikka_id:tä",id1,"ei löydy paikkataulusta!"))
  }
  if (!(id2 %in% kaikki_paikat$paikka_id)) {
    warning(paste("Paikka_id:tä",id2,"ei löydy paikkataulusta!"))
  }
  if (all(!(c(id1,id2) %in% kaikki_paikat$paikka_id))) {
    warning("Putkipari ohitettu")
  } else {
    pai1 <- kaikki_paikat[kaikki_paikat$paikka_id == id1, ]
    pai2 <- kaikki_paikat[kaikki_paikat$paikka_id == id2, ]
    
    if ((pai1$lon != pai2$lon) | (pai1$lat != pai2$lat))
      cat("DB ero:",pai1$tunnus,pai1$paikka_id,"¤", pai2$tunnus,pai2$paikka_id,"\n")
    if ((!pai1$yhd_lon %in% pai2$yhd_lon) | !(pai1$yhd_lat %in% pai2$yhd_lat))
      cat("ELYkrd ero:",pai1$tunnus,pai1$paikka_id,"¤",pai2$tunnus,pai2$paikka_id,"\n")
  }
}
rm(id1,id2, pai1,pai2)

```


```{r save-elevation-table}
#### Save the collated pipe top and ground elevation file
# Perform last cheks
elevtbl_sanitycheck(putki_yhd_temp,
                    skip_yla_tyyppi = c("Routaputki","Lumikeppi","Lysimetri"))
elevtbl_sanitycheck(putki_yhd_trim,
                    skip_yla_tyyppi = c("Routaputki","Lumikeppi","Lysimetri"))

# Select only required cols
korkeudet_outfile <-
  putki_yhd_trim[ ,c("paikka_id","tunnus","asema_tunnus","asema_nimi",
                     "PaikkaTyyppi_Id","yla","yla_source","yla_verif",
                     "maa","maa_source","maa_verif","KM2_laatu")]
#Round numeric cols to two decimal digits
korkeudet_outfile <-
  dplyr::mutate_if(korkeudet_outfile, is.numeric, ~round(., 2)) |>
  dplyr::mutate(PaikkaTyyppi_Id = as.numeric(PaikkaTyyppi_Id))
# Order into a sensible order
korkeudet_outfile <-
  korkeudet_outfile[stringr::str_order(korkeudet_outfile$tunnus, numeric = T), ]
korkeudet_outfile <- korkeudet_outfile[order(korkeudet_outfile$asema_tunnus), ]
# Warn of IDs which would be uploaded but are not included in the investigated tbl
checked <- data.table::fread(f.path(D$data,"Miralle_DRAFT_korkeudet_CHECKED.csv"))$paikka_id
for (id in korkeudet_outfile$paikka_id) {
  if (!(id %in% checked)) warning("ID ",id," puuttuu tarkistettujen taulusta!")
}; rm(checked)

# Write the Muonio 1305 station test case
data.table::fwrite(
  korkeudet_outfile[korkeudet_outfile$asema_tunnus %in% "1305", ],
  check_dir(D$tulokset,"korkeudet","korkeudet_tietokantaan_MUONIO.csv"),
  dec = ".", sep = ";", bom = T)

# Write the results to a file (Sans MUONIO)
data.table::fwrite(
  korkeudet_outfile[!(korkeudet_outfile$asema_tunnus %in% "1305"), ],
  check_dir(D$tulokset,"korkeudet","korkeudet_tietokantaan_DRAFT.csv"),
  dec = ".", sep = ";", bom = T)


# Write also the full resulting dataset. 
# Includes values for places, which were not wanted in the large db update
# (also extra info: ground frost pipes top elevs, etc. )
putki_yhd_temp <- putki_yhd_temp[stringr::str_order(putki_yhd_temp$tunnus, numeric=T), ]
putki_yhd_temp <- putki_yhd_temp[order(putki_yhd_temp$asema_tunnus), ]
data.table::fwrite(putki_yhd_temp,
  check_dir(D$tulokset,"korkeudet","korkeudet_FULL_DATA.csv"),
  dec = ".", sep = ";", bom = T)

# Write also a table which includes all alternative ground / top elevs
all_kork <-
  merge(putki_yhd_temp,
        kaikki_paikat[,c("paikka_id","yhd_yla","yhd_maa","dtm_elev","km2_elev")],
        by = "paikka_id")
all_kork["dtm_yla"] <- all_kork$dtm_elev + all_kork$putki 
all_kork["km2_yla"] <- all_kork$km2_elev + all_kork$putki

all_kork <- all_kork[stringr::str_order(all_kork$tunnus, numeric = T), ]
all_kork <- all_kork[order(all_kork$asema_tunnus), ]
data.table::fwrite(all_kork,
  check_dir(D$tulokset,"korkeudet","korkeudet_FULLEST_DATA.csv"),
  dec = ".", sep = ";", bom = T)


# Write the dropped places in a file (both db upload compliant and full infos)
putki_yhd_dropped <- putki_yhd_dropped[stringr::str_order(putki_yhd_dropped$tunnus, numeric=T), ]
putki_yhd_dropped <- putki_yhd_dropped[order(putki_yhd_dropped$asema_tunnus), ]
korkeus_dropped_out <-
  putki_yhd_dropped[ ,c("paikka_id","tunnus","asema_tunnus","asema_nimi",
                     "PaikkaTyyppi_Id","yla","yla_source","yla_verif",
                     "maa","maa_source","maa_verif","KM2_laatu")]
korkeus_dropped_out <- korkeus_dropped_out[!korkeus_dropped_out$PaikkaTyyppi_Id %in% "Lumikeppi",]
data.table::fwrite(korkeus_dropped_out,
  check_dir(D$tulokset,"korkeudet","korkeudet_tietokantaan_DROPPED.csv"),
  dec = ".", sep = ";", bom = T)

all_dropped <-
  merge(putki_yhd_dropped,
        kaikki_paikat[,c("paikka_id","yhd_yla","yhd_maa","dtm_elev","km2_elev")],
        by = "paikka_id")
all_dropped["dtm_yla"] <- all_dropped$dtm_elev + all_dropped$putki 
all_dropped["km2_yla"] <- all_dropped$km2_elev + all_dropped$putki
all_dropped <- all_dropped[stringr::str_order(all_dropped$tunnus, numeric = T), ]
all_dropped <- all_dropped[order(all_dropped$asema_tunnus), ]
data.table::fwrite(all_dropped,
  check_dir(D$tulokset,"korkeudet","korkeudet_DROPPED_FULL_DATA.csv"),
  dec = ".", sep = ";", bom = T)


```

``` {r plot-new-elevs}

# Create a table for plotting, calculate required values,
# replacing missing db values with 0 in a separate col
putki_yhd_plot <- putki_yhd_temp
putki_yhd_plot[["maa_delta"]] <- putki_yhd_plot$maa - putki_yhd_plot$KorkeusMaa
putki_yhd_plot[["maa_filled"]] <- putki_yhd_plot$maa -
  replace(putki_yhd_plot$KorkeusMaa, is.na(putki_yhd_plot$KorkeusMaa), 0)
putki_yhd_plot[["maa_ver"]] <-
  ifelse(putki_yhd_plot$maa-putki_yhd_plot$KorkeusMaa==0, 0, NA)
putki_yhd_plot[["yla_delta"]] <- putki_yhd_plot$yla - putki_yhd_plot$yla_db
putki_yhd_plot[["yla_filled"]] <- putki_yhd_plot$yla -
  replace(putki_yhd_plot$yla_db, is.na(putki_yhd_plot$yla_db), 0)
putki_yhd_plot[["yla_ver"]] <-
  ifelse(putki_yhd_plot$yla-putki_yhd_plot$yla_db==0, 0, NA)

# Drop cases from plots which need manual review
putki_yhd_plot <- putki_yhd_plot[
  !(putki_yhd_plot$maa_verif %in% c("CHECKS_FAILED") |
    putki_yhd_plot$yla_verif %in%
      c("CHECKS_FAILED","BAD_DB?(NO_VERIF)","ELY_IN_DB?",p0("ELY±",l_ely,"m"))), ]

# Plot the overall picture of the differences in elevations
maa_deltaplot <- ggplot(putki_yhd_plot) +
  geom_histogram(aes(maa_filled, fill=alpha("black",.7)), binwidth=2) +
  geom_histogram(aes(maa_delta, fill=alpha("red",.7)), na.rm=T, binwidth=2) +
  geom_histogram(aes(maa_ver, fill=alpha("yellow",.9)), na.rm=T, binwidth=2) +
  coord_trans(y="pseudo_log") +
  scale_y_continuous(breaks=c(0,5,10,25,50,100,200,500,1000)) +
  theme_classic() + theme(legend.position="top") +
  labs(x="Ground elevation change [m]", y="Count") +
  scale_fill_manual(
    name=NULL,
    values=alpha(c("black","red","yellow"), c(.7,.7,.9)),
    labels=c("Originally missing elevs","Existing elevs",
             "Unchanged")) +
  guides(fill = guide_legend(reverse=T))

yla_deltaplot <- ggplot(putki_yhd_plot) +
  geom_histogram(aes(yla_filled, fill=alpha("black",.7)), binwidth=2) +
  geom_histogram(aes(yla_delta, fill=alpha("red",0.7)), na.rm=T, binwidth=2) +
  geom_histogram(aes(yla_ver, fill=alpha("yellow",.9)), na.rm=T, binwidth=2) +
  coord_trans(y="pseudo_log") +
  scale_y_continuous(breaks=c(0,5,10,25,50,100,200,500,1000)) +
  theme_classic() + theme(legend.position="top") +
  labs(x="Top elevation change [m]", y="Count") +
  scale_fill_manual(
    name=NULL,
    values=alpha(c("black","red","yellow"), c(.7,.7,.9)),
    labels=c("Originally missing elevs","Existing elevs",
             "Unchanged")) +
  guides(fill = guide_legend(reverse=T))

save2(maa_deltaplot, f.path(D$tulokset,"korkeudet","maanpinta_muutos"))
save2(yla_deltaplot, f.path(D$tulokset,"korkeudet","ylapaa_muutos"))


# Get the results into tables for presenting
top_elev_sources <- data.frame(
  source = c(
    ref_meas = sum(putki_yhd_plot$yla_source %in% c("Pietari")),
    db = sum(putki_yhd_plot$yla_source %in% c("DB")),
    dtm = sum(putki_yhd_plot$yla_source %in% c("DTM+DB","DTM+ELY","DTM+MISSING")),
    km2 = sum(putki_yhd_plot$yla_source %in% c("KM2+DB","KM2+ELY","KM2+MISSING")),
    ely = sum(putki_yhd_plot$yla_source %in% c("ELY")))
)
top_elev_sources[["frac"]] <-
  top_elev_sources$source / sum(!is.na(putki_yhd_plot$yla_source)) * 100
# No. of places where structure top elevation value changed:
sum(top_elev_sources$source) - top_elev_sources["db","source"]
.Last.value / sum(top_elev_sources$source)


ground_elev_sources <- data.frame(
  source = c(
    ref_meas = sum(putki_yhd_plot$maa_source %in%
                     c("Pietari","Pietari-Pietari", "DB-Pietari")),
    db = sum(putki_yhd_plot$maa_source %in% c("DB","DB-MISSING")),
    dtm = sum(putki_yhd_plot$maa_source %in%
                c("DTM","DTM+DB-DB","DTM+MISSING-MISSING","DTM+ELY-ELY")),
    km2 = sum(putki_yhd_plot$maa_source %in%
                c("KM2","KM2+ELY-ELY","KM2+DB-DB","KM2+MISSING-MISSING")),
    ely = sum(putki_yhd_plot$maa_source %in% c("ELY","DB-ELY")))
)
ground_elev_sources[["frac"]] <-
  ground_elev_sources$source / sum(!is.na(putki_yhd_plot$maa_source)) * 100
# No. of places where ground elevation value changed:
sum(ground_elev_sources$source) - ground_elev_sources["db","source"]
.Last.value / sum(ground_elev_sources$source)

# Count how many values were changed from DB value, and how many DB values were verified
cat("Average change of changed ground values:",
mean(replace(putki_yhd_temp$maa-putki_yhd_temp$KorkeusMaa,
             putki_yhd_temp$maa-putki_yhd_temp$KorkeusMaa == 0, NA), na.rm=T),
"\nAverage change of changed structure top values:",
mean(replace(putki_yhd_temp$yla-putki_yhd_temp$yla_db,
             putki_yhd_temp$yla-putki_yhd_temp$yla_db == 0, NA), na.rm=T)
)


#rm(putki_yhd_temp)

```



```{r session-info, results='hide'}

sessinfo_path <- f.path("doc","R_session_info.txt") #Define session info outfile

ops <- options()
options(width=80L) # Ensure that the printout formatting remains constant

sessioninfo::session_info(to_file=sessinfo_path) # Output the session info
# Take the resulting file and omit the library path
modif_info <- readLines(sessinfo_path) |> gsub(".*C:/Users/.*", "", x=_)
cat(modif_info, file=sessinfo_path, sep="\n"); rm(modif_info) #Write the modified file

# Append also base::sessionInfo()
cat("\n+++ RESULTS OF base::sessionInfo:", capture.output(sessionInfo()),"",
    file=sessinfo_path, sep="\n", append=T) #("fill=N" is inadequate for formatting a list)

options(ops); rm(ops) # Restore options
rm(sessinfo_path)

```