#install.packages('rgdal')
#install.packages('dplyr')
#install.packages('sf')

require(rgdal)
require(dplyr)
library(sf)

################################################################################################
################################################################################################
#
# Title       : SOM2Map.R
# Author      : Brittany M. Krzyzanowski
#
################################################################################################
################################################################################################


#Input MatLab SOM table output named regiontable1.csv
umat=read.csv("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\regiontable1.csv")
str(umat)
utm15nCRS=st_crs(poly)

#plot lat long csv to shp

umat_pt <- st_as_sf(umat, coords = c("x", "y"), crs = utm15nCRS)
plot(umat_pt$geometry,
    main = "Map of Plot Locations")

st_write(umat_pt,"C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk\\umat_pt.shp", driver = "ESRI Shapefile")
U<- readOGR("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "umat_pt") 
  

#theissen polygon function
voronoipolygons <- function(x) {
  require(deldir)
  require(sp)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  z <- deldir(crds[,1], crds[,2])
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
    y=crds[,2], row.names=sapply(slot(SP, 'polygons'), 
    function(x) slot(x, 'ID'))))
}


v <- voronoipolygons(U)
plot(v)

writeOGR(v,"C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "v_poly",driver="ESRI Shapefile",overwrite_layer=TRUE)  
vpoly<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "v_poly") 

#umat<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "umat_pt") 

proj4string(U) <- CRS("+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83")
writeOGR(U,"C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "umat",driver="ESRI Shapefile",overwrite_layer=TRUE)  
umat<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "umat") 

v<- readOGR("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "v_poly")  
proj4string(v) <- CRS("+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83")
writeOGR(v,"C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "v_poly2",driver="ESRI Shapefile",overwrite_layer=TRUE)  
v<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "v_poly2") 

join <- st_join(v, left = FALSE, umat["umat"])
plot(join)

st_write(join, "C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk\\join_regions.shp",update=TRUE)
v_join<- readOGR("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "join_regions") 
v_join<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "join_regions")  

poly<- readOGR("C:\\Users\\krzyzanowski\\Desktop\\Step2_R", "metro_bg_poly") 
plot(poly)

poly<- readOGR("C:\\Users\\krzyzanowski\\Desktop\\Step2_R", "metro_bg_poly")  
proj4string(poly) <- CRS("+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83")
writeOGR(poly,"C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "poly2",driver="ESRI Shapefile",overwrite_layer=TRUE)  
poly<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "poly2") 

m=st_join(poly, v_join, join = st_nearest_feature)
st_write(m, "C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk\\m.shp",update=TRUE)
m<- st_read("C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\Junk", "m") 
 
#plot(m)
m$x <- NULL
m$y <- NULL
m$count <- NULL
m$GeoX <- NULL
m$GeoY <- NULL
m$TOTPOP= as.numeric(as.character(m$TOTPOP_CY))

st_write(m, "C:\\Users\\krzyzanowski\\Desktop\\Step2_R\\output\\umat_map.shp",update=TRUE)

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