#path_to_coordinates <- "/Users/Beniamino/Downloads/NFDB_point_txt/"
path_to_data <- "/Users/Beniamino/Downloads/NFDB_point_txt/NFDB_point_20150911_large_fires.txt"
library(data.table)

# read original dataset
datasettxt <- data.table(read.table(path_to_data, sep=",", header = TRUE))
# select only interesting columns: when, where, how much
data <- datasettxt[,c(5,6,7,8,9,14), with=FALSE]
# select only years common to every region (>=1980)
data <- data[YEAR_>=1980] 
# remap to 0.5 degrees precision-> close fire happened at the same place
for (i in 1:length(data[,LATITUDE])) {
  whole <- floor(data[i,LATITUDE])
  decimal <- (data[i,LATITUDE]-floor(data[i,LATITUDE]))
  if (decimal>0.25 & decimal <= 0.75 ) {data[i,"LATITUDE"]<-(whole+0.75)}
  else {
    if (decimal<0.5) {data[i,"LATITUDE"]<-(whole+0.25)}
    else {data[i,"LATITUDE"]<-(whole+1.25)}
  }
  
  whole <- floor(data[i,LONGITUDE])
  decimal <- (data[i,LONGITUDE]-floor(data[i,LONGITUDE]))
  if (decimal>0.25 & decimal <= 0.75 ) {data[i,"LONGITUDE"]<-(whole+0.25)}
  else {
    if (decimal<0.5) {data[i,"LONGITUDE"]<-(whole-0.25)}
    else {data[i,"LONGITUDE"]<-(whole+0.75)}
  }
}

# create list with data in 1980-2014 counting number of occurences of fire at same location, for every year and for evey latlon location
lens <- list()
for (i in 1980:2014) {
len <- list()
idx <- 1
  for (j in 1:53) {
   for (k in 1:144) {
    len[[idx]] <- nrow(data[YEAR_==i & LATITUDE ==(43.25+0.5*j) & LONGITUDE == (-69.25-0.5*k)]) # the number of rows corresponds to number of fires
    idx <- idx +1
  }
  }
rl <- rbind(len)
lens[[i]] <- t(rl)
}

# add list containing the locations in the same order
lon <- list()
lat <- list()
idx <- 1
for (j in 1:53) {
  for (k in 1:144) {
    lon[[idx]] <- (-69.25-0.5*k)
    lat[[idx]] <- (43.25+0.5*j)
    idx <- idx +1
  }
}
lat <- t(rbind(lat))
lon <- t(rbind(lon))
# create new data.table containing the occurences of fire per location (rows) per year (columns)
len_table <- data.table(t(rbindlist(lens)))
for (i in 1980:2014) {
  setnames(len_table,paste0("V",(i-1979)),paste0("YEAR_",i))
  }
# add locations 
len_table[,LAT:=lat]
len_table[,LON:=lon]

# add the sum of all fires per location in the timeinterval and the average number per location
sumlist <- list()
avglist <- list()
tot <- 0
for (i in 1:nrow(len_table)) {
  sums <- sum(len_table[i,1:35,with=FALSE])
  sumlist[[i]] <- sums
  avglist[[i]] <- sums/35
  tot <- tot + sums
}

sumcol <- t(rbind(sumlist))
len_table[,SUM:=sumcol]
len_table[, AVG:=t(rbind(avglist))]

# select only points north of 45 degrees
len_table <- len_table[LAT>=45.]

# let's try to plot it as it is
# load libraries
library(stringr)
library(ggmap)
# Get Map
map_Location <- get_map(location = c(lon = mean(as.numeric(len_table$LON)), 
                                     lat = mean(as.numeric(len_table$LAT))) ,
                        zoom = 3, language = "en-EN",
                        source = "google", maptype = "terrain")

# Plot all map points using the GDD0 for the hue
tempPlot <- ggmap(map_Location) + geom_point(data = len_table,
                                             aes(x = as.numeric(LON), y = as.numeric(LAT), color=as.numeric(SUM)), size=2) + 
  xlab("Longitude") + ylab("Latitude") + ggtitle("Number of Fires over Canada")
print(tempPlot) 



canadian_wildfires <- data.table(as.numeric(len_table[,LAT]))
canadian_wildfires[,LON:=as.numeric(len_table$LON)]
canadian_wildfires[,SUM:=as.numeric(len_table$SUM)]
canadian_wildfires[,AVG:=as.numeric(len_table$AVG)]
setnames(canadian_wildfires,"V1","LAT")

path_to_save <- "/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/Dataframes/"
write.csv(canadian_wildfires,file=paste0(path_to_save,"canadian_wildfires_1980-2014.csv"), quote = FALSE)


dataset_na_east <- read.csv("/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/boreal_na_east_database.txt",sep = " " , header = FALSE)
setnames(dataset_na_east,"V11","LON")
setnames(dataset_na_east,"V12","LAT")
setnames(dataset_na_east,"V6","FF")



counter <- 0
counter_missing <- 0 
for (i in 1:length(dataset_na_east[,"FF"])) {
  if (nrow(canadian_wildfires[LAT==dataset_na_east[i,"LAT"] & LON==dataset_na_east[i,"LON"]]) == 1) {counter <- counter +1}
  if (nrow(canadian_wildfires[LAT==dataset_na_east[i,"LAT"] & LON==dataset_na_east[i,"LON"]]) == 0) {counter_missing <- counter_missing +1}
  
}

for (i in 1:length(dataset_na_east[,"FF"])) {
  if (nrow(canadian_wildfires[LAT==dataset_na_east[i,"LAT"] & LON==dataset_na_east[i,"LON"]]) == 1) {
    maxvalue <- max(dataset_na_east[i,"FF"], canadian_wildfires[LAT==dataset_na_east[i,"LAT"] & LON==dataset_na_east[i,"LON"],AVG] )
    dataset_na_east[i,"FF"] <- maxvalue}
}

write.csv(dataset_na_east,file = "/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/boreal_na_east_database_canadianupdate.txt",sep = " ", col.names = FALSE, row.names = FALSE)

dataset_na_west <- read.csv("/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/boreal_na_west_database.txt",sep = " " , header = FALSE)
setnames(dataset_na_west,"V11","LON")
setnames(dataset_na_west,"V12","LAT")
setnames(dataset_na_west,"V6","FF")



counter <- 0
counter_missing <- 0 
for (i in 1:length(dataset_na_west[,"FF"])) {
  if (nrow(canadian_wildfires[LAT==dataset_na_west[i,"LAT"] & LON==dataset_na_west[i,"LON"]]) == 1) {counter <- counter +1}
  if (nrow(canadian_wildfires[LAT==dataset_na_west[i,"LAT"] & LON==dataset_na_west[i,"LON"]]) == 0) {counter_missing <- counter_missing +1}
  
}

for (i in 1:length(dataset_na_west[,"FF"])) {
  if (nrow(canadian_wildfires[LAT==dataset_na_west[i,"LAT"] & LON==dataset_na_west[i,"LON"]]) == 1) {
    maxvalue <- max(dataset_na_west[i,"FF"], canadian_wildfires[LAT==dataset_na_west[i,"LAT"] & LON==dataset_na_west[i,"LON"],AVG] )
    dataset_na_west[i,"FF"] <- maxvalue}
}

write.csv(dataset_na_west,file = "/Users/Beniamino/Documents/OneDrive/MPI-M/iPython_Notebooks/Boreal_Area/boreal_na_west_database_canadianupdate.txt",sep = " ", col.names = FALSE, row.names = FALSE)






for (i in 1901:1902) {
  thaw_dpt <- data.table(read.table(paste0(path_to_data, i,".txt")))
  idx <-i-1900
  mean_thaw_dpt<- data.table(thaw_dpt$V2)
  setnames(mean_thaw_dpt,paste0("V1"), paste0("MEAN_THAW_DEPTH_",i))
  
  setnames(thaw_dpt,paste0("V1"), paste0("MEDIAN_THAW_DEPTH_",i))
  setnames(thaw_dpt,paste0("V2"), paste0("MEAN_THAW_DEPTH_",i))
  setnames(thaw_dpt,paste0("V3"), paste0("STD_DEV_THAW_DEPTH_",i))
  setnames(thaw_dpt,paste0("V4"), paste0("MIN_THAW_DEPTH_",i))
  setnames(thaw_dpt,paste0("V5"), paste0("MAX_THAW_DEPTH_",i))
  
  
  
  mean_thaw_dpt[,LON:=locs_global$LON]
  mean_thaw_dpt[,LAT:=locs_global$LAT]
  setkey(mean_thaw_dpt,LON,LAT)
  
  
  #locs_global <- locs_global[thaw_dpt]
  #thaw_list[[idx]] <- thaw_dpt
  locs_global <- locs_global[mean_thaw_dpt]
  print(i)
}





write.csv(locs_global,file=paste0(path_to_data,"mean_1901-2002.csv"), quote = FALSE)
#write.table(locs_global,file=paste0(path_to_data,"1901-2002.txt"),sep = "\t")
#locs_global_reordered <- data.table(read.table(paste0(path_to_data, "locs_global_reordered.csv"), sep=","))
#setnames(locs_global_reordered,c("V1", "V2", "V3"), c("LAT", "LON", "THAW"))

library(proj4)
data <- cbind(locs_global_reordered$LON, locs_global_reordered$LAT, locs_global_reordered$THAW)
#src <- +proj=laea  +lat_0=90.0 +lon_0=0.0 +x_0=0 +y_0=0 +a=6371228.0 +units=m
#dest <- +init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
unproj <- ptransform(data, '+proj=laea  +lat_0=90.0 +lon_0=0.0 +x_0=0 +y_0=0 +units=km','+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0', silent=TRUE)
library(ggplot2)
library(ggmap)
unproj.table <- as.data.frame(unproj)
map_Location <- get_map(location = c(lon = mean(unproj.table$V2), 
                                     lat = mean(unproj.table$V1)) ,
                        zoom = 1, language = "en-EN",
                        source = "google", maptype = "terrain")
tempPlot <- ggmap(map_Location) + geom_point(data = unproj.table,
                                             aes(x = V2, y = V1, colour=V3), size=5) + 
  xlab("Longitude") + ylab("Latitude") + ggtitle("Test")
print(tempPlot)