list.of.packages <- c("gdata","R.utils","rgdal","raster","rbioclim","stringr","stringi",
                      "proj4","data.table","FSA","plyr","wordcloud","KernSmooth","gdtools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages) > 0) install.packages(new.packages)

lapply(list.of.packages, FUN = function(X) {
    do.call("require", list(X)) 
})

devtools::install_github("MoisesExpositoAlonso/rbioclim") # To download WorldClim
library(rbioclim) # To download WorldClim

path <- getwd()
setwd(path)
set.seed(1601720)
rasterOptions(maxmemory=1.1e9)

# Set up necessary directories
if (!dir.exists("WorldClim")) {dir.create("WorldClim")}
if (!dir.exists("DEM")) {dir.create("DEM")}

source("Translator.R")

clim_p <- recursive.getData("pres",path="WorldClim",var="bio")
clim_p <- clim_p$pres[[1]]
plot(clim_p)

clim_f <- recursive.getData("CC2650",path="WorldClim",var="bio")
clim_f <- clim_f$CC2650[[1]]
plot(clim_f)

if (file.exists("Provided_DEM.tif")) {
    dem <- raster("Provided_DEM.tif")
} else {
    if (!file.exists("DEM/15-C.tif")) {
  download.file("http://www.viewfinderpanoramas.org/DEM/TIF15/15-C.zip", 
                destfile="DEM/isl_2010.zip")
  unzip("DEM/isl_2010.zip",exdir="DEM")
        }
dem <- raster("DEM/15-C.tif",sep="")
    }
plot(dem)

writeRaster(dem,"Iceland_DEM.tif",type="GTiff",overwrite=TRUE)
writeRaster(clim_f,"Iceland_Clim_Future.tif",type="GTiff",overwrite=TRUE)
writeRaster(clim_p,"Iceland_Clim_Present.tif",type="GTiff",overwrite=TRUE)

dem <- raster("Iceland_DEM.tif")
clim_f <- raster("Iceland_Clim_Future.tif")
clim_p <- raster("Iceland_Clim_Present.tif")

unlink("topo",recursive=TRUE)
unlink("WorldClim",recursive=TRUE)

#proj <- CRS("+init=epsg:9040 +datum=WGS84")
proj <- CRS("+proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000
                    +datum=WGS84 +units=m +no_defs")

ext <- raster()  # Dummy target extent raster
extent(ext) <- extent(2400000,3000000,4300000,4700000)
crs(ext) <- proj

ext_proj <- projectExtent(ext,dem) # The DEM
dem <- crop(dem,ext_proj)
dem[dem <= 0] <- NA # Since the ocean is reported as having Z = 0---we want it as NA

ext_proj <- projectExtent(ext,clim_f) # The WorldClim data. 
clim_f <- crop(clim_f,ext_proj)
clim_p <- crop(clim_p,ext_proj)

if (!file.exists("Provided_DEM.tif")) {
    dem <- projectRaster(dem,crs=crs(ext),res=200)
    dem <- crop(dem,ext)
}
plot(dem)

clim_f <- projectRaster(clim_f,crs=crs(ext),res=1957) #Project
clim_p <- projectRaster(clim_p,crs=crs(ext),res=1957)

clim_f <- resample(clim_f,dem) # Resample
clim_p <- resample(clim_p,dem)

clim_f <- crop(clim_f,ext) # and Crop
clim_p <- crop(clim_p,ext)

plot(clim_f)

writeRaster(dem,"Iceland_DEM.tif",type="GTiff",overwrite=TRUE)
writeRaster(clim_f,"Iceland_Clim_Future.tif",type="GTiff",overwrite=TRUE)
writeRaster(clim_p,"Iceland_Clim_Present.tif",type="GTiff",overwrite=TRUE)


if (file.exists("Provided_Results.gz")) {
    arc <- fread("Provided_Results.gz", encoding="UTF-8")
} else {
    arc <- arc.dataconvert("results.json.gz")
}

setnames(arc,names(arc),c("NUMBER","ID","Y","X","DATASET","CATEGORY","ENTRY","CONCEPT","COMBINATOR","RELATED","CONTEXT"))
arc[, names(arc) := lapply(.SD,trimws,which='left')]

tail(arc)

arc[,c("X","Y") := lapply(.SD,as.numeric), .SDcols = c("X","Y")]

arc[, `:=`(X = X * ..pi / 180, Y = Y * ..pi / 180)]

# Now we can convert to our projection system
arc[,`:=`(X = ptransform(.(X,Y),src.proj=CRS("+proj=longlat +datum=WGS84"),
                         dst.proj=proj)[[1]],
          Y = ptransform(.(X,Y),src.proj=CRS("+proj=longlat +datum=WGS84"),
                         dst.proj=proj)[[2]])]
arc[1:25,.(ID,X,Y)]

output <- arc[,lapply(.SD,mean), by=ID, .SDcols = c("X","Y")][,.(X,Y)]
output <- na.omit(output)

output <- SpatialPointsDataFrame(output[,.(X,Y)],output,proj4string=proj)
plot(dem)
plot(output,add=TRUE)

#mapview::mapview(output)

concepts <- read.csv("CONCEPT_DATA_FRAME.csv",stringsAsFactors = FALSE)
concepts <- as.data.table(concepts)

query <- concepts[concepts$NAME == "butchery",HASH] 

output <- arc[str_detect(CONCEPT,query), # Filters by hash, gets XY
              lapply(.SD,mean), by=.(ID,CONCEPT), .SDcols = c("X","Y")][,.(X,Y,CONCEPT)]
output <- na.omit(output) # Drops missing values
output <- SpatialPoints(output[,.(X,Y)],proj4string=proj) # Adds projection information
#mapview::mapview(output, map.types = "OpenTopoMap")

plot(dem)
plot(output, add=TRUE)

query <- c("agricultur| farm")

hash <- concepts[str_detect(concepts$NAME, query)] # Finds the right hashes
hashes <- paste(hash$HASH,collapse="|") # Puts them in a format readable by stringr
output <- arc[str_detect(CONCEPT,hashes), # Filters by hash, gets XY and concept
              lapply(.SD,mean), by=.(ID,CONCEPT), .SDcols = c("X","Y")][,.(X,Y,CONCEPT)]

output <- output[, .(CONCEPT = unlist(tstrsplit(CONCEPT,":"))), by=c("X","Y")] # Convert to List
output <- output[CONCEPT %in% hash$HASH,.(CONCEPT),by=c("X","Y")] # Filter unwanted concepts

output[, `:=`(DUMMY = match(CONCEPT,..hash$HASH))]
output <- na.omit(output) # Drops missing values
output <- SpatialPointsDataFrame(output[,.(X,Y)],output[,.(CONCEPT,DUMMY)],
                                 proj4string=proj) # Adds projection information
output$CONCEPT <- hash$NAME[match(output$CONCEPT, hash$HASH)]

# with MapView
#mapview::mapview(output, zcols = "CONCEPT", map.types = "OpenTopoMap") # Plot the result

# with plot
# Plot the DEM again to clear the display, then plot the points
plot(dem)
plot(output, pch=(1:length(output))[output@data$DUMMY],
     col=(1:length(output))[output@data$DUMMY],
     add= TRUE)

# Add a legend
legend("topright",hash$NAME,
       pch=(1:length(output)),
       col=(1:length(output)))

plot(clim_p)

clim_p <- clim_p / 10

output$Temp <- extract(clim_p,output)

hist(output$Temp~output$CONCEPT, 
     same.ylim=FALSE, w=2.5,
     xlab = "Temperature (degrees Celcius)",
     ylab = "Frequency")

clim_f <- clim_f / 10 # Don't forget to convert the predicted temperatures to degrees!
change <- clim_f - clim_p
plot(change)

output <- arc[,lapply(.SD,mean), by=.(ID,CONCEPT), .SDcols = c("X","Y")][,.(X,Y,CONCEPT)] # Combine Sites
output <- na.omit(output) # Exclude null values
output <- SpatialPointsDataFrame(output[,.(X,Y)],output[,.(X,Y,CONCEPT)],proj4string=proj) # Convert to spatial
output$TempChange <- extract(change,output) # Extract the temperature change
warming_sites <- output[complete.cases(output@data$TempChange),] # Keep only sites with real values
warming_sites <- warming_sites[warming_sites$TempChange >= 2.0,] # Only keep sites warming by more than 2 deg

output$Elevation <- extract(dem,output) # Extract elevation
flooding_sites <- output[complete.cases(output@data$Elevation),] # Keep only sites with real values
flooding_sites <- flooding_sites[flooding_sites$Elevation <= 1,] # Only keep sites below 1 meter above sea level

datasets <- read.csv("DATASET_HASH_OUT.csv",
                     stringsAsFactors = F,
                     encoding='UTF-8') # Get dataset hash list
datasets <- as.data.table(datasets)
datasets[, names(datasets) := lapply(.SD,trimws,which='left')]
query <- "Threat" # Set our keyword
threat <- datasets[str_detect(NAME,query),HASH] # Pull out the correct hash
threat <- toupper(threat) # Since arc hashes are in all caps
threatened_sites <- arc[DATASET == threat, 
                        lapply(.SD,mean), by=.(ID,CONCEPT),
                        .SDcols = c("X","Y")][,.(X,Y,CONCEPT)]
threatened_sites <- SpatialPointsDataFrame(threatened_sites[,.(X,Y)],
                                          threatened_sites,
                                          proj4string = proj)

cols = c("X","Y","CONCEPT","TYPE") # Define the columns we'll later keep
warming_sites$TYPE <- "Temperature" # Assign risk types
flooding_sites$TYPE <- "Flooding"
threatened_sites$TYPE <- "DataARC"
threatened <- data.table(rbind(warming_sites@data[,(cols)], # Combine risk types
                               flooding_sites@data[,(cols)],
                               threatened_sites@data[,(cols)]))
threatened[, c("X","Y") := lapply(.SD,as.numeric),.SDcols = c("X","Y")]
threatened_points <- SpatialPointsDataFrame(threatened[,.(X,Y)],threatened, proj4string = proj) # Convert to Spatial
#mapview::mapview(threatened_points,zcol="TYPE") # Plot

plot(dem)
plot(warming_sites, pch=1, col=1,add=T)
plot(flooding_sites, pch=2, col=2,add=T)
plot(threatened_sites, pch=3, col=3,add=T)
legend("topright", c("Temperature Change","Flooding","DataARC 'At-Risk'"),
       pch=1:3,col=1:3)

threatened <- threatened[, .(CONCEPT = unlist(tstrsplit(CONCEPT,":"))), 
                         by=c("X","Y","TYPE")]
threatened[, CONCEPT := ..concepts$NAME[match(CONCEPT,..concepts$HASH)]]
threatened <- na.omit(threatened)

counts <- count(threatened$CONCEPT)
wordcloud(counts$x,counts$freq,random.order = FALSE, colors=brewer.pal(8, "Dark2"))