library(raster)
library(cartography)
library(sf)
library(SpatialPosition)

mtq <- st_read(system.file("shape/martinique.shp", package="cartography"))
# use WGS84 proj
mtq_latlon <- st_transform(mtq, 4326)
# this call throw an error but seems to work...
getData('SRTM', lon=-61, lat=14)
# import raster
ras <- raster("srtm_24_10.tif")
# crop on martinique area
mtq_ras <- crop(ras, st_bbox(mtq_latlon)[c(1,3,2,4)])
# aggregate the raster
mtq_ras <- aggregate(mtq_ras, fact=4,fun=mean)
mtq_ras <- projectRaster(mtq_ras,  crs = "+proj=utm +zone=20 +datum=WGS84 +units=m +no_defs")
# break values
bv <- c(seq(0,1300,100),1339)
# contour extraction
mtq_cont <- rasterToContourPoly(r = mtq_ras, breaks = bv, mask = as(mtq, "Spatial"))
# custom palette
pal <- c("#5D9D52", "#8DBC80", "#B8D9A9", "#FDEBBE", "#F7E0AC", "#F2D69B", 
         "#EDCC8A", "#E8C279", "#E2B563", "#DBA84C", "#D49B36", "#BA8428", 
         "#9A6A1E", "#7B5114")
# sp to sf
k <- st_as_sf(mtq_cont)
# order the sf
k <- k[order(k$center),]
png(filename = "mtq.png", width = 700, height = 800, res = 100, bg = NA)
par(mar = c(0,0,1.2,0))
plot(st_geometry(mtq), col = NA, border = NA, bg = "lightblue")
for(i in 1:nrow(k)){
  p <- st_geometry(k[i,])
  plot(p + c(-50, 50), add=T, border = "#ffffff90",col = "#ffffff90")
  plot(p + c(100, -100),  col = "#00000090", add=T, border = "#00000090")
  plot(p, col = pal[i], border = "NA", add=T)  
}
legendChoro(pos = "left", breaks = bv, col = pal, nodata = F,
            title.txt = "Elevation\n(metres)", cex = 0.8)
layoutLayer(title = "Martinique Relief", north = T,
            sources = 'T. Giraud, 2018', author = "SRTM, 2018", col = "lightblue", 
            tabtitle = T, coltitle = "black")
dev.off()