# Load required libraries
setwd("V:/UHI/models/R_2022_training")
library(matrixStats)
library(quantreg)
library(RobustLinearReg)

#setwd("C:/Users/kpostma/OneDrive - DOI/NLCDShrub/20240226_Hua_UHI_R")

# Read CSV file
# data <- read.csv("V:/UHI/analysis/census_analys/data/csv_input/sea_bg_data_2020.csv")
data <- read.csv("V:/UHI/CGI/analysis/wuh/input/wuh_cgi_tmean_y_20_23.csv")
# data <- read.csv("V:/UHI/CGI/analysis/wuh/input/wuh_cgi_tmin_y_20_23.csv")
# data <- read.csv("V:/UHI/CGI/analysis/wuh/input/wuh_cgi_tmax_y_20_23.csv")

names(data)

# Extract x and y variables from data
# xs <- c('md_n_2020211','vp_n_2020211', 'md_d_2020211','vp_d_2020211', 'ls_0729_2020')
xs <- c('CGI_LST_2020','MODIS_2020')
# xs <- c('CGI_LST_2023','MODIS_2023')
# ys <- c('md_d_2020211','vp_d_2020211', 'ec_n_2020212')
ys <- c('MODIS_2020','VIIRS_2020')
# ys <- c('MODIS_2023','VIIRS_2023')

run_model <- function(xname, yname) {
  #xname<-'MeanLST'
  #yname<-'HOUSE'
  #xname<-'INTS_MeanLST'
  #yname<-'HOUSE'
  x <- as.double(unlist(data[xname]))
  y <- as.double(unlist(data[yname]))
  final <- data.frame(x = x, y = y)
  final <- final[complete.cases(final), ]

  # Fit linear regression to get slope, R square, and P-value
  #lm_fit <- lm(y ~ x)
  #lm_intercept <- coef(lm_fit)[1]
  #lm_slope <- coef(lm_fit)[2]
  #lm_r_squared <- summary(lm_fit)$r.squared
  #lm_p_value <- summary(lm_fit)$coefficients[2, 4]

  fit <- theil_sen_regression(formula = y ~ x, data = final)
  ts_intercept <- coef(fit)[1]
  ts_slope <- coef(fit)[2]
  #pearson_r <- cor(x = final$y, y = predict(fit, final))
  #r_squared <- pearson_r * pearson_r
  ts_r_squared <- summary(fit)$r.squared
  ts_p_value <- summary(fit)$coefficients[2, 4]

  # Plot the data and the regression line
  plot(x, y, main = "Scatter plot with Theil-Sen Estimation",
       xlab = paste(xname), ylab = paste(yname))
  
  #abline(coef = c(lm_intercept, lm_slope), col = 'blue', lwd = 2)
  abline(coef = c(ts_intercept, ts_slope), col = 'red', lwd = 2)
  #legend("topleft", legend = c(paste0('LM Slope: ', round(lm_slope, 4), '\nR square: ', round(lm_r_squared, 4), '\nP-value: ', round(lm_p_value, 4), '\n'),
                               #paste0('Theil-Sen Slope: ', round(ts_slope, 4), '\nR square: ', round(ts_r_squared, 4), '\nP-value: ', round(ts_p_value, 4), '\n')),
         #col = c('blue','green'), lty = 1, cex = 0.8)
  dev.copy(png, paste0('plot','_X_',xname,'_Y_',yname,'.png'), units = 'px', width = 2048, height = 2048, res = 300)
  dev.off()
}

#run_model('MeanLST','HOUSE')

for (xname in xs) {
  for (yname in ys) {
    run_model(xname,yname)
  }
}

