mountain_parameter.R

pacman::p_load_gh("uribo/fgdr", "franapoli/pbarETA")
pacman::p_load(sf, raster, tabularaster, tidyverse, broom, ggthemes, plotly, Metrics)

plot_surface <- function(mat) plot_ly(x=mat[, 1], y=colnames(mat), z=mat[, -1]) %>% add_surface()

options(pillar.sigfig = 10)
list_xml <- list.files("data", pattern="FG-GML-5133-36.+DEM5.+.xml", full.names=T)
pb <- txtProgressBar(min=2, max=length(list_xml), style=3, char="█", width=getOption("width") - 30)
source <- list()
for(i in 1:length(list_xml)){
  source[[i]] <- read_fgd_dem(list_xml[i], resolution=5, return_class="raster") %>%  setValues(., getValues(.) %>% if_else(. == -9999, -1, .)) %>%
    as_tibble(xy=T) %>% mutate(source=list_xml[i])
  setTxtProgressBar(pb, i)
}

df <- bind_rows(source) %>% rename(lat=y, lon=x)
df <- df %>% filter(dplyr::between(lat, 34.265, 34.285)) %>% filter(dplyr::between(lon, 133.835, 133.857))
df %>% summary
df <- df %>% group_by(lat, lon) %>% summarise(alt=mean(cellvalue, na.rm=T)) %>% ungroup
df %>% summary

df <- df %>% mutate(z=scale(alt) %>% as.numeric, x=scale(lon) %>% as.numeric, y=scale(lat) %>% as.numeric)

ggplot(df, aes(x=lon, y=lat, z=alt)) + geom_contour(bins=20, aes(color=stat(level))) +
  coord_equal() + labs(x="longitude", y="latitude", color="altitude") + theme_pander() + labs(title="Contour of Mt. Iino")



##### Plateau distribution ####
# plot univariate plateau dist.
dplateau <- function(x, loc=0, scale=1, alpha=1){
  return(alpha/pi * sin(pi/(2*alpha))/(1 + ((x - loc)/scale)^(2*alpha)))
}

pmap(tibble(a=1:6),
     function(a) {stat_function(data=. %>% mutate(a=a, label=as.character(a)), aes(color=label), fun=dplateau, args=list(alpha=a))}
) %>% 
  reduce(.init=ggplot(data=tibble(x=seq(-4, 4, .1)), aes(x=x)), `+`) +
  labs(color=expression(alpha), y="density", title="Density of Standard Plateau Dist.") +
  scale_color_pander() + theme_pander()

pmap(tibble(mean=seq(-1, 1, .4)) %>% mutate(label=factor(mean)),
     function(mean, label) {stat_function(data=. %>% mutate(mean=mean, label=label), aes(color=label), fun=dplateau, args=list(loc=mean))}
) %>% 
  reduce(.init=ggplot(data=tibble(x=seq(-4, 4, .1)), aes(x=x)), `+`) +
  labs(color="location\nparameter", y="density", title="Density of Standard Plateau Dist.") + theme_pander() +
  scale_color_pander(breaks=seq(-1, 1, .4))

pmap(tibble(s=seq(.5, 2, .4)),
     function(s) {stat_function(data=. %>% mutate(s=s), aes(color=as.character(s)), fun=dplateau, args=list(scale=s))}
) %>% 
  reduce(.init=ggplot(data=tibble(x=seq(-4, 4, .1)), aes(x=x)), `+`) +
  labs(color="scale\nparameter", y="density", title="Density of Standard Plateau Dist.") +
  scale_color_pander() + theme_pander()

dmvplateau <- function(x, y, m_x, m_y, s_x, s_y, a){
  a^2 * sin(pi/(2*a))^2 / (pi^2 * s_x * s_y) * ((1 + ((x - m_x)/s_x)^(2*a)) * (1 + ((y-m_y)/s_y)^(2*a)))^-1
}
# plotly で 3D プロット
plot_surface <- function(mat) plot_ly(x=mat[, 1], y=colnames(mat), z=mat[, -1]) %>% add_surface()

expand.grid(seq(-4, 4, .1), seq(-4, 4, .1)) %>% as_tibble %>% set_names(c("x", "y")) %>% mutate(density = dmvplateau(x, y, 0, 0, 2, 2, 2)) %>%
  spread(key=y, value=density) %>% as.matrix %>% plot_surface %>%
  layout(title="Density of 2D Plateau Distribution", scene=list(zaxis=list(title="density")))

##### fit parametric curves ####

param_std <- df %>% dplyr::select(lat, lon) %>% summarise(mean_lat=mean(lat), sd_lat=sd(lat), mean_lon=mean(lon), sd_lon=sd(lon))

eq_2d_Gaussian <- z ~ k*(2*pi*s_x*s_y*(1-rho^2)^(1/2))^(-1) *
  exp(-(2 - 2 * rho^2)^-1 * (((x - m_x)/s_x)^2 + ((y - m_y)/s_y)^2 - 2*rho*(x - m_x)*(y - m_y)/(s_x * s_y)))
#eq_2d_Student <- z ~ k * gamma(v/2+1)/(gamma(v/2) * v * pi * s_x * s_y * (1-rho^2)^(1/2)) *
#  ( 1 + (v * (2 - 2 * rho^2)) * ((x - m_x)^2/s_x^2 + (y - m_y)^2/s_y^2 - 2*rho*(x - m_x)*(y - m_y)/(s_x * s_y)))^-(v/2 + 1)
eq_2d_Cauchy <- z ~ k * gamma(3/2)/(gamma(1/2) * pi * s_x * s_y * (1-rho^2)^(1/2)) *
  ( 1 + (2 - 2 * rho^2) * ((x - m_x)^2/s_x^2 + (y - m_y)^2/s_y^2 - 2*rho*(x - m_x)*(y - m_y)/(s_x * s_y)))^-(3/2)
eq_2d_Student5 <- z ~ k * gamma(7/2)/(gamma(5/2) * 5 * pi * s_x * s_y * (1-rho^2)^(1/2)) *
  ( 1 + (5 * (2 - 2 * rho^2)) * ((x - m_x)^2/s_x^2 + (y - m_y)^2/s_y^2 - 2*rho*(x - m_x)*(y - m_y)/(s_x * s_y)))^-(7/2)
eq_2d_logis <- z ~ k * 2 * exp(- (x - m_x)/s_x - (y - m_x)/s_y) / (1 + exp(-(x - m_x)/s_x) + exp(-(y - m_y)/s_y))^3
# eq_2d_plateau <- z ~ k * a^2 * sin(pi/(2*a))^2 / ((pi^2 * s_x * s_y) * ((1 + ((x - m_x)/s_x)^(2*a)) * (1 + ((y - m_y)/s_y)^(2*a))))
eq_2d_plateau1 <- z ~ k * 2 * sin(pi/2)^2 / ((pi^2 * s_x * s_y) * ((1 + ((x - m_x)/s_x)^2) * (1 + ((y - m_y)/s_y)^2)))
eq_2d_plateau2 <- z ~ k * 4 * sin(pi/4)^2 / ((pi^2 * s_x * s_y) * ((1 + ((x - m_x)/s_x)^4) * (1 + ((y - m_y)/s_y)^4)))

# test 
nls(formula = eq_2d_plateau2, data=df,
    start=list(m_x=0, m_y=0, s_x=sd(df$x), s_y=sd(df$y), k=25),
    lower=c(m_x=-100, m_y=-100, s_x=0, s_y=0, k=0),
    algorithm="port",
    trace=T, control=nls.control(maxiter=100, minFactor=1/10000, warnOnly=T))

inits <- list(m_x=0, m_y=0, s_x=sd(df$x), s_y=sd(df$y), k=15)
lower <- c(m_x=-100, m_y=-100, s_x=-.0001, s_y=0, k=0)
edit <-function(li, name, val){
  li[[name]] <- val
  return(li)
}
equations <- c(Gaussian=eq_2d_Gaussian, Cauchy=eq_2d_Cauchy, `Student (5)`=eq_2d_Student5,
               Logistic=eq_2d_logis, `Plateau (1)`=eq_2d_plateau1, `Plateau (2)`=eq_2d_plateau2)
set.seed(42)
result <- pmap(
  list(
    EQ=equations,
    INITS=list(c(inits, rho=0), c(edit(inits, "k", 30), rho=0), c(inits, rho=0), inits, edit(inits, "k", 25), edit(inits, "k", 25))
    ),
  function(EQ, INITS) nls(EQ, data=df, start=INITS,
                          lower=c(m_x=-100, m_y=-100, s_x=0, s_y=0, k=0),
                          algorithm="port",
                          control=nls.control(maxiter=100, warnOnly=T, minFactor=1/10000))) %>%
  set_names(names(equations))

#### Post Diagnosis ####

# original estimates
(df_result <- map2(result, names(result), function(RESULT, NAME) tidy(RESULT) %>% mutate(name=NAME)) %>%
  bind_rows %>% dplyr::select(name, term, estimate) %>% spread(key=term, value=estimate))

(df_metrics <- tibble(
  name=names(result),
  RMSE=map_dbl(result, function(x) sqrt(mean(resid(x)^2))),
  MAE=map_dbl(result, function(x) mean(abs(resid(x))) )
) %>% arrange(RMSE)) # %>% xtable::xtable(digits=4) %>% print(include.rownames=F)

df_result %>% arrange(name) %>%
  mutate(m_x=param_std$mean_lon + param_std$sd_lon*m_x, m_y=param_std$mean_lat + param_std$sd_lat*m_y,
         s_x=abs(s_x)*param_std$sd_lon, s_y=abs(s_y)*param_std$sd_lat) # %>% xtable::xtable(digits=4) %>% print(include.rownames=F)

# residual plot

map2(result, names(result), function(OBJ, NAME){df %>% mutate(resid= resid(OBJ)) %>% arrange(x, y) %>%
  ggplot(aes(x=x, y=y, fill=resid)) + geom_tile() + theme_pander() + labs(title=paste("Residuals of ", NAME))}
)

# plotly で 3D プロット
plot_surface <- function(mat) plot_ly(x=mat[, 1], y=colnames(mat), z=mat[, -1]) %>% add_surface()

map2(result, names(result), function(OBJ, NAME) df %>% dplyr::select(x, y) %>% mutate(resid=resid(OBJ)) %>% spread(key=y, value=resid) %>% as.matrix() %>% plot_surface %>% layout(title=paste("Residuals of", NAME), scene=list(zaxis=list(title="residual"))))