Computing Linear Shape Metric

per.area_atlas.R

library(rgdal)
library(spatstat)
library(maptools)
library(rgeos)

setwd("D:/METHODES/1211_ELONGATEDCOUNTRIES/")

countries <- readOGR("D:/DATAS/NATURAL EARTH/ne_110m_admin_0_countries_lakes/ne_110m_admin_0_countries_lakes.shp", "ne_110m_admin_0_countries_lakes")


##############
# ELONGATION #
##############

fragstats.circle <- function(x) {
els <- sapply(1:nrow(x), function(i) {
  ow <- as(countries[i, ], "owin")
  inc <- incircle(ow)
  dsc <- disc(inc$r, c(inc$x, inc$y))
  elgt <- gArea(as(dsc, "SpatialPolygons")) / gArea(countries[i, ])
  return(elgt) 
})
return(els)
}


################
# CALCULATIONS #
################

elgtds <- fragstats.circle(countries)


##########
# COLORS #
##########

cols <- list('Sub-Saharan Africa' = rgb(253, 180, 98, maxColorValue=255),
             'Middle East & North Africa' = rgb(255, 255, 179, maxColorValue=255),
             'Antarctica' = rgb(141, 211, 199, maxColorValue=255),
             'South Asia' = rgb(190, 186, 218, maxColorValue=255), 
             'Europe & Central Asia' = rgb(179, 222, 105, maxColorValue=255), 
             'North America' = rgb(128, 177, 211, maxColorValue=255),
             'Latin America & Caribbean' = rgb(251, 128, 114, maxColorValue=255),
             'East Asia & Pacific' = rgb(252, 205, 229, maxColorValue=255))


###########
# MATRICE #
###########

# configuration de la matrice
nrows <- round(sqrt(nrow(countries)))
ncols <- ifelse(nrows^2 < nrow(countries), nrows + 1, nrows)

# couleurs des pays
colors <- sapply(countries$region_wb, function(x) cols[[as.character(x)]])

# création de l'image
png(file="IMG/elongation.png", width=3000, height=3000)
par(mfrow=c(nrows,ncols), bg=gray(.8), oma=c(1,1,1,1))
sapply(rev(order(elgtds)), function(x) {
  plot(countries[x, ], type="n")
  plot(countries, col=gray(.7), border=gray(.6), add=T)
  plot(countries[x, ], col=colors[[x]], border=NA, add=T) # pays
  title(strwrap(iconv(countries$name[x], "UTF-8", "latin1"), width=15, prefix="\n", initial=""), col.main=colors[[x]], cex.main=3.5) # titre
})
dev.off()


##################
# CARTE DU MONDE #
##################

png(file="IMG/monde.png", width=3000, height=1500)
par(bg=gray(.8))
plot(countries, col=colors, border=gray(.7), lwd=1)
legend("bottomleft", 
       names(cols), cex=3, 
       fill=as.character(cols), border=NA,
       bty="o", box.col=gray(.7));
dev.off()