ptompalski

library(lidR)
library(tidyverse)
library(glue)
library(ggforce) #geom_circle
library(gifski)

#set your paths here
dir_out <- "C:/Users/ptompals/OneDrive - NRCan RNCan/__WORK/graphics"

#temp dir will be a subfolder of the main dir

ptompalski

# function to plot a crossection of a point cloud
# works with lidR package

plot_crossection <- function(las,
                             p1 = c(min(las@data$X), mean(las@data$Y)),
                             p2 = c(max(las@data$X), mean(las@data$Y)),
                             width = 4, colour_by = NULL)
{
  colour_by <- enquo(colour_by)
  data_clip <- clip_transect(las, p1, p2, width)

ptompalski

library(raster)
library(tidyverse)
library(lubridate)
library(insol)
library(av)

#input data is a digital surface model (DSM) generated with airborne lidar
dsm1 <- raster("dsm.tif")

#convert to martix (required by rayshader)

ptompalski

bcf <- function(observed, model_name) {
  # another way to correct for back-transformation bias
  # based on this:
  # https://stats.stackexchange.com/questions/361618/how-to-back-transform-a-log-transformed-regression-model-in-r-with-bias-correcti
  
  # 1. Compute exp(Xβ^), i.e. the retransformed but unadjusted prediction
  a <- exp(predict(model_name))
  
  # 2. Regress Y against exp(Xβ^) without an intercept. Call the resulting regression coefficient γ.
  b <- lm(observed ~ a -1) #adding -1 in the formula removes the intercept

ptompalski

# command line based raster mosaicing

# example usage:
# mosaic.py c:\rasters\to\mosaic\*.asc c:\output\raster\mosaic.tif

# script requires ArcMap license


import arcpy
import os

ptompalski

mode <- function(x) {
  if(is.numeric(x)) {  x <- round(x,2)}
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}