library(readr)
library(ggplot2)
library(scales)
library(dplyr)
library(sf)
library(rgdal)
library(data.table)
library(terra)
#####################
#data

velky <- read_csv("D:/data_GEDI/export_clip_Z/velky/finale_lite.csv")


polygons1km <- vect('D:/data_GEDI/mriz/M1km.shp')
polygons500m <- vect('D:/data_GEDI/mriz/M500m.shp')
polygons250m <- vect('D:/data_GEDI/mriz/M250m.shp')
#settings

neighbors <- c(9)
thresholds <- c(5)

# filtr nastav  & num_detectedmodes <= 4 
model <- velky %>%
  filter((quality_flag == TRUE & rx_assess_flag == 0 & num_detectedmodes <= 4))

####################

polygons1km <- project(polygons1km, "EPSG:4326")
polygons500m <- project(polygons500m, "EPSG:4326")
polygons250m <- project(polygons250m, "EPSG:4326")
model <- st_as_sf(model, coords = c("XX", "YY"), crs = 4326)
coords <- st_coordinates(model)
results_df <- data.frame(neighbor = integer(),
                         threshold = integer(),
                         RMSE = numeric(),
                         count = integer())

for (n in neighbors) {
  
  # Find 10 nearest neighbors (including the point itself, so k = 11 for 10 others)
  nn <- RANN::nn2(coords, coords, k = n)
  
  # Calculate median height for the neighbors
  model$median_height_neighbors <- apply(nn$nn.idx, 1, function(idx) {
    median(model$elev_lowestmode[idx])
  })
  
  for (t in thresholds) {
    
    
    model1 <- model %>%
      filter(abs(elev_lowestmode - median_height_neighbors) <= t)
    
    squared_differences <- model1$diff^2
    mean_squared_error <- mean(squared_differences)
    rmse <- sqrt(mean_squared_error)
    
    count <- nrow(model1)
    results_df <- rbind(results_df, data.frame(neighbor = n, threshold = t, RMSE = rmse, count = count))
    print(t)
    print(n)
    print(Sys.time())
    print(rmse)
  }
}
model1 <- model1[model1$diff < 30, ]
coords <- st_coordinates(model1)
model1 <- data.frame(model1, X = coords[,1], Y = coords[,2])
model1$geometry <- NULL
model1 <- vect(model1, geom=c("X", "Y"), crs="EPSG:4326")
##250m
modelint <- intersect(model1, polygons250m)
modelint <- as.data.frame(modelint, geom="XY")
point_counts <- modelint %>%
  group_by(Id) %>%
  summarise(count = n())
polygons <- as.data.frame(polygons250m)
joined <- full_join(polygons, point_counts, by = "Id")

print(count(modelint))
print(summary(joined))
print("250m")
##500m
modelint <- intersect(model1, polygons500m)
modelint <- as.data.frame(modelint, geom="XY")
point_counts <- modelint %>%
  group_by(Id) %>%
  summarise(count = n())
polygons <- as.data.frame(polygons500m)
joined <- full_join(polygons, point_counts, by = "Id")

print(count(modelint))
print(summary(joined))
print("500m")
##1km
modelint <- intersect(model1, polygons1km)
modelint <- as.data.frame(modelint, geom="XY")
point_counts <- modelint %>%
  group_by(Id) %>%
  summarise(count = n())
polygons <- as.data.frame(polygons1km)
joined <- full_join(polygons, point_counts, by = "Id")

print(count(modelint))
print(summary(joined))
print("1km")