library(rgdal)
library(rgeos)
library(raster)
library(ggplot2)
library(ggthemes)
library(ggsn)
library(ggmap)
library(maps)
library(maptools)
library(grid)
library(GISTools)
library(celestial)
library(bibtex)
library(plyr) # this needs to be always before dplyr, otherwise, it's masking important functions
library(dplyr)
library(SoDA)

par(bg = "white")
puerto.ayora = data.frame(lon = -90.3138000,
                          lat = -0.7401800,
                          row.names = "puerto.ayora")
el.garrapatero = data.frame(lon = -90.221744,
                            lat = -0.687647,
                            row.names = "el.garrapatero")
locations = rbind(puerto.ayora,
                  el.garrapatero)
coord.deg = locations
class(coord.deg)
coordinates(coord.deg)<-~lon+lat
class(coord.deg)
proj4string(coord.deg) # nope
proj4string(coord.deg)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
URL <- "https://osm2.cartodb.com/api/v2/sql?filename=public.galapagos_islands&q=select+*+from+public.galapagos_islands&format=geojson&bounds=&api_key="
fil <- "gal.json"
if (!file.exists(fil)) download.file(URL,
                                     fil)
gal <- readOGR(fil,
               "OGRGeoJSON")

URL.topo <- "https://osm2.carto.com/api/v2/sql?filename=contour_1&q=select+*+from+public.contour_1&format=geojson&bounds=&api_key="
fil.topo = "public.geojson.json.topo"
if (!file.exists(fil.topo)) download.file(URL.topo,
                                          destfile = fil.topo)
topo <- readOGR("public.geojson.json.topo",
                layer ="OGRGeoJSON")
topo <- spTransform(topo,
                    proj4string(gal)) # epsg:4326 is decimal degree,
                                      # epsg:31983 is meters only,
                                      # 31969 works too
a <- aggregate(gal)
b <- disaggregate(a)
# Setting the theme
theme_set(theme_bw())
#Create a custom color scale
coord.deg<-spTransform(coord.deg,
                       CRS(proj4string(gal)))
# double check that they match
identical(proj4string(coord.deg),
          proj4string(gal)) # look also proj4string(topo)
id=1:length(coord.deg)
my_pts <- SpatialPointsDataFrame(coords = coord.deg,
                                 data=data.frame(id=1:length(coord.deg)))

gal@data$id <- seq(gal@data$fid) - 1
gal_union <- unionSpatialPolygons(gal,
                                  IDs = rep(1,length(gal)))
gal_map <- fortify(gal)
gal_ff <- fortify(gal)
gal_union_ff <- fortify(gal_union)

topo_map <- fortify(topo)
topo_map <- plyr::join(topo_map,
                 topo@data, by="id")
names(topo_map)
# ggplot can't deal with a SpatialPointsDataFrame so we
# can convert back to a data.frame
my_pts <- data.frame(my_pts)
my_pts.final = my_pts[ ,2:3]
# we're not dealing with lat/long but with x/y
# this is not necessary but for clarity change variable names
names(my_pts.final)[names(my_pts.final)=="lat"]<-"y"
names(my_pts.final)[names(my_pts.final)=="lon"]<-"x"

# Complete Galápagos Island -------------------------------------------

### settings ###

##  colour of the map
# ground
ground = "darkseagreen" # "#8FBC8F"
ground.contour = "black"
ground.contour.size = .75 # around the islands

# water
water = "lightsteelblue1" # "#CAE1FF" # you can put it white "#FFFFFF" if you want
# water = NA # suggestion: NA => no color

# border contour
border.map = "black"
panel.border.size = 1
# GPS points
gps.points.col = c("#FF3030") # "red"
size.of.points = 0.7
# grid colour
inside.grid.c = "blue"
inside.grid.s = 0.1
inside.grid.t = "dotted"

# scale colour for topographic maps
topo.low = "#8FBC8F" # "darkseagreen"
topo.high = "darkgreen" # "#006400"

#scale bar size & location
dist.in.km = 50 # half of the distance of the bar
location.scale = "bottomleft"

# legend position
leg.pos = "right" # if "none" removes all legends

## text
# Title and axis
title.graph = "Îles Galápagos"
y.lab.map = "Latitude"
x.lab.map = "Longitude"
legend.title = "Altitude"
ticks.length = .2
# size (note that the scalebar is different)
text.s = 12
title.mult = 1.5 # multiplicator relative to text.s, so that the title is bigger
# scale
text.scale = text.s/3

## topographic information
# lines
topo.line.size = 0.1
transparent.topo = 0.4
linetype.topo = 1
### end settings ###

gg1 <- ggplot()
gg <- gg1 + geom_path(aes(x = long,
                          y = lat,
                          group = group),
                      data = gal_map,
                      colour = ground.contour, # will draw each islands
                      size = ground.contour.size) + # Size!
  geom_polygon(data = gal_map,
               aes(x = long,
                   y = lat,
                   group = group),
               color = ground, # colour of the ground
               fill = ground) # colour of the space between the polygon. Otherwise, you get weird lines.

gg <- gg + ggsn:::scalebar(data = gal_map, # scale bar
                           dist= dist.in.km, # set the size of 1/2 of the scalebar
                           location= location.scale,
                           model = "WGS84", # options ("WGS84", "GRS80", "Airy", "International", "Clarke", "GRS67")
                           dd2km = TRUE, # this must be true in order to use the distance as km
                           st.size=text.scale, # number to indicate the scale bar text size. It is passed to the size argument of annotate function.
                           st.bottom = TRUE) # scale bar text is displayed at the bottom of the scale bar
gg <- gg + #theme_map() + # veryyyyyyyyyyyyyyyyyy sober
  theme(plot.title = element_text(family="Helvetica",
                                  size = text.s*title.mult,
                                  face = "bold"),
        legend.position = leg.pos,
        legend.text = element_text(size = text.s,
                                   colour = "black",
                                   angle = 0),
        axis.text = element_text(family="Helvetica",
                                 size=text.s), # text size
        axis.title = element_text(family="Helvetica",
                                  size=text.s), # text size
        axis.ticks.length = unit(ticks.length,
                                 "cm"), # length of the axis ticks
        panel.background = element_rect(fill = water), # white for water
                # remove the default grey, you can add blue color for
                # water (# background making illusion that this is water)
        panel.grid.major = element_blank(), #element_line(colour = inside.grid.c,
                                            # linetype = inside.grid.t,size = inside.grid.s),
                                            # element_blank() if you don't want the grid
        legend.background = element_blank(), # remove background from legend
        panel.grid.major = element_blank(), # remove grid
        panel.grid.minor = element_blank(), # remove grid
        strip.background = element_blank(), # remove grid
        plot.background = element_blank(), # remove white background of
            # the plot (if you want to put your graph in a publication or presentation that has another colour)
        panel.border = element_rect(colour = border.map, # border to the map
                                    fill=NA,
                                    size=panel.border.size))
gg <- gg +  #geom_polygon(data = subset.c, aes(x = long, y = lat, fill = subset.c[,"id"], group = group)) +
  coord_equal()

gg <- gg + geom_path(data = topo_map, # adding the elevation line
                     aes(long,
                         lat,
                         group=group, # group by island
                         colour = elev),
                     size = topo.line.size , # based on your feeling !
                     alpha=transparent.topo, # if you want to see through the lines (depth perseption)
                     linetype = linetype.topo) + # 1 - solid line
  scale_color_gradient(low = topo.low, # colour lower elevation
                       high = topo.high, # colour higher elevation
                       guide = "colourbar") + # guide = "none" if you don't want
  labs(title = title.graph, # title of the graph
       y = y.lab.map, # title of y axis
       x = x.lab.map, # title of x axis
       colour = legend.title) #+
# Adding points only at the end
gg <- gg + geom_point(data=my_pts.final[,], # sites points
                      aes(x = x,
                          y = y,
                          size = size.of.points),
                      color = gps.points.col) +
  geom_text(data=my_pts.final[1,], aes(x = x, y = y,
                                       label=c("Puerto Ayora")),
            size = text.scale,
            vjust = +0.6, # moves things down with a +
            hjust = -0.1) +
  geom_text(data=my_pts.final[2,], aes(x = x, y = y,
                                       label=c("El Garrapatero")),
            size = text.scale, vjust = +0.2, hjust = -0.1)


gg <- gg + coord_equal() # to get a nice projection
df <- data.frame(x1 = -89.5, x2 = -89.5, y1 = 1, y2 = 1.25)