#!/usr/bin/Rscript
## gadm_plot.r -- create SVG plot of administrative areas of a country
## usage: ./gadm_plot.r <primary> [<secondary>+]
## where <primary> and <secondary> are the three-letter country codes as known to GADM
## see http://gadm.org/
##
## = Other arguments =
## -numbers : display region ID numbers at centre point
## -labels : display surrounding country labels at centre point
## -level2 : display Level 2 boundaries (if available)
## -resolution <float> : point displacement resolution (degrees, or equivalent units)
## -minpoints <int> : only show polygons with more than <int> points
##
## Author: David Eccles (gringer) 2010 <programming@gringer.org>
setupPlot <- function(polydata, bounds = NA, waterColour = "#C6ECFF"){
require(sp);
if(is.na(bounds)){
bounds <- polydata@bbox;
}
countryName <- as.character(polydata@data$NAME_FAO[1]);
cat("Setting up plot boundaries ...", file = stderr());
bbox.points <- rbind(as.vector(bounds)[c(1,2)],
as.vector(bounds)[c(1,4)],
as.vector(bounds)[c(3,4)],
as.vector(bounds)[c(3,2)],
as.vector(bounds)[c(1,2)]);
par(mar = c(0,0,0,0), bg = waterColour);
## adjust aspect ratio to fix Y distortion at mean latitude (i.e. a quick mercator-ish projection)
## a more correct approach would be to convert the coordinates to a particular projection
plot(bbox.points,
type = "n", asp = 1/cos(mean(bounds["r2",]) * pi/180),
bty = "n", axes = FALSE, xlab = "", ylab = "",
xaxs = "r", yaxs = "r");
cat("done!\n", file = stderr());
}
level0plot <- function(polydata, lineWidth = 2,
landColour = "#FEFEE9", waterColour = "#C6ECFF", resolution = 0.001, minPoints = 0){
require(sp);
countryName <- as.character(polydata@data$NAME_FAO[1]);
cat("Drawing Level 0 for",countryName,"...", file = stderr());
for(p in polydata@plotOrder){
polysubplot(polydata@polygons[p][[1]], regionColour = landColour, holeColour = waterColour, lineWidth = lineWidth,
resolution = resolution, minPoints = minPoints);
}
cat("done!\n", file = stderr());
}
level1plot <- function(polydata, colours = NULL, nameDisplacement = 0.2, landColour = "#FEFEE9",
waterColour = "#C6ECFF", drawNames = TRUE, drawNumbers = FALSE, resolution = 0.001, minPoints = 0){
countryName <- as.character(polydata@data$NAME_0[1]);
cat("Drawing Level 1 for",countryName,"...", file = stderr());
require(sp);
for(p in polydata@plotOrder){
if((!is.null(colours)) && (length(colours) == length(polydata@plotOrder))){
landColour = colours[p];
}
if(polydata@data$ENGTYPE_1[p] == "Water body"){
polysubplot(polydata@polygons[p][[1]], regionColour = waterColour, holeColour = landColour,
resolution = resolution, minPoints = minPoints);
} else {
polysubplot(polydata@polygons[p][[1]], regionColour = landColour, holeColour = waterColour,
resolution = resolution, minPoints = minPoints);
}
}
if(drawNames || drawNumbers){
## conversions for unicode characters [please add others as necessary]
polydata@data$NAME_1 <- sub("\xfd","\u00fd",polydata@data$NAME_1); # y + ?acute
polydata@data$NAME_1 <- sub("\xfc","\u00fc",polydata@data$NAME_1); # u + diaresis
polydata@data$NAME_1 <- sub("\xe2","\u00e2",polydata@data$NAME_1); # a + circumflex
polydata@data$NAME_1 <- sub("\xed","\u00ed",polydata@data$NAME_1); # ?
polydata@data$NAME_1 <- sub("\xe9","\u00e9",polydata@data$NAME_1); # e + acute
polydata@data$NAME_1 <- sub("\xe8","\u00e8",polydata@data$NAME_1); # e + grave
polydata@data$NAME_1 <- sub("\xe1","\u00e1",polydata@data$NAME_1); # ?
polydata@data$NAME_1 <- sub("\xda","\u00da",polydata@data$NAME_1); # ?
polydata@data$NAME_1 <- sub("\xc1","\u00c1",polydata@data$NAME_1); # A + ?acute
range <- apply(bbox(polydata),1,diff);
for(p in sort(polydata@plotOrder)){
centrePoint <- polydata@polygons[p][[1]]@labpt;
plotPoint <- centrePoint;
nd <- nameDisplacement;
if(nd > 0){
plotPoint[1] <- plotPoint[1] + runif(1) * nd * range[1] - range[1] * nd / 2;
plotPoint[2] <- plotPoint[2] + runif(1) * nd * range[2] - range[2] * nd / 2;
lines(rbind(centrePoint, plotPoint));
}
if(drawNames){
text(x = t(plotPoint), labels = polydata@data$NAME_1[p]);
}
if(drawNumbers){
cat(": '''",as.character(p), "''' (", polydata@data$NAME_1[p], ")\n", sep="", file = stdout());
text(x = t(plotPoint), labels = as.character(p));
}
}
}
cat("done!\n", file = stderr());
}
level2plot <- function(polydata, colours = NULL, bounds = polydata@bbox, resolution = 0.001){
countryName <- as.character(polydata@data$NAME_0[1]);
cat("Drawing Level 2 for",countryName,"...", file = stderr());
require(sp);
landColour = NA;
waterColour = NA;
for(p in polydata@plotOrder){
if((!is.null(colours)) && (length(colours) == length(polydata@plotOrder))){
landColour = colours[p];
}
polysubplot(polydata@polygons[p][[1]], lineWidth = 0.5, regionColour = landColour, holeColour = waterColour,
resolution = resolution);
}
cat("done!\n", file = stderr());
}
level3plot <- function(polydata, colours = NULL, bounds = polydata@bbox, resolution = 0.001){
countryName <- as.character(polydata@data$NAME_0[1]);
cat("Drawing Level 2 for",countryName,"...", file = stderr());
require(sp);
landColour = NA;
waterColour = NA;
for(p in polydata@plotOrder){
if((!is.null(colours)) && (length(colours) == length(polydata@plotOrder))){
landColour = colours[p];
}
polysubplot(polydata@polygons[p][[1]], lineWidth = 0.25, regionColour = landColour, holeColour = waterColour,
resolution = resolution);
}
cat("done!\n", file = stderr());
}
linedist <- function(point.matrix){
## calculates minimum distance between point (xP,yP) and line [(x1,y1) - (x2,y2)]
## assumes input is as follows
## x1 y1
## xP yP
## x2 y2
if(!all(dim(point.matrix) - c(3,2) == 0)){
stop("input must be a 3x2 matrix");
}
## distance formula from http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
x1 <- point.matrix[1,1];
y1 <- point.matrix[1,2];
xP <- point.matrix[2,1];
yP <- point.matrix[2,2];
x2 <- point.matrix[3,1];
y2 <- point.matrix[3,2];
if((x1 == x2) && (y1 == y2)){
## exclude case where denominator is zero
return(sqrt((x1-xP)^2 + (y1 - yP)^2));
}
dist = abs((x2 - x1)*(y1 - yP) - (x1 - xP)*(y2 - y1)) /
sqrt((x2 - x1)^2 + (y2 - y1)^2);
return(dist);
}
angle.point.extract <- function(points.matrix){
mdim <- dim(points.matrix);
if(mdim[2] != 2){
stop("must be an mx2 matrix");
}
if(mdim[1] < 3){
stop("must be at least 3x2 matrix");
}
return(aperm(array(rbind(t(points.matrix[1:(mdim[1] - 2),]),
t(points.matrix[2:(mdim[1] - 1),]),
t(points.matrix[3:(mdim[1] - 0),])),
c(2,3,(mdim[1] - 2))),c(2,1,3)));
}
polysubplot <- function(polydata, regionColour = "#FEFEE9", holeColour = "#C6ECFF",
lineColour = "#646464", lineWidth = 1, resolution = 0.001,
minPoints = 0){
require(sp);
for(p in polydata@plotOrder){
coords <- polydata@Polygons[p][[1]]@coords;
## transform into desired projection
## remove points below resolution limit
pointCount <- dim(coords)[1];
if(pointCount > 2){
diffs <- apply(angle.point.extract(coords),3,linedist);
keeps <- rep(FALSE,dim(coords)[1]);
keeps[1] <- TRUE;
change <- 0;
for(x in 1:length(diffs)){
change <- change + diffs[x];
if(change > resolution){
keeps[x+1] <- TRUE;
change <- 0;
}
}
keeps[pointCount] = TRUE;
coords <- coords[keeps,];
}
## only two points, separated by less than the resolution
if((dim(coords)[1] == 2) && (sqrt(sum(diff(coords)^2)) < resolution)){
coords <- coords[1,];
}
pointCount <- dim(coords)[1];
if((is.null(dim(coords))) || (dim(coords)[1] == 1)){
if(minPoints <= 1){
if(is.null(dim(coords))){
points(t(coords), pch = 21, cex = resolution / 0.1 * lineWidth, col = lineColour, bg = regionColour);
} else {
points(coords, pch = 21, cex = resolution / 0.1 * lineWidth, col = lineColour, bg = regionColour);
}
}
} else {
## now draw image
if(pointCount >= minPoints){
if(polydata@Polygons[p][[1]]@hole){
polygon(coords, col = holeColour, lwd = lineWidth, border = lineColour);
} else {
polygon(coords, col = regionColour, lwd = lineWidth, border = lineColour);
}
}
}
}
}
drawPrimary0 <- function(countryData, colours = NULL, resolution = 0.001,
landColour = "#FEFEE9", waterColour = "#C6ECFF", drawNames = FALSE, minPoints = 0){
require(sp);
if(is.character(countryData)){
cat("Loading Level 0 data for", countryCode, "(primary country)...", file = stderr());
con <- url(paste("http://gadm.org/data/rda/",countryCode,"_adm0.RData",sep=""));
load(con);
close(con);
cat("done!\n", file = stderr());
countryBounds <- bbox(gadm);
level0plot(gadm, landColour = landColour, waterColour = waterColour, lineWidth = 2,
resolution = resolution, minPoints = minPoints);
} else {
gadm <- countryData;
countryBounds <- bbox(gadm);
level0plot(gadm, landColour = landColour, waterColour = waterColour, lineWidth = 2,
resolution = resolution, minPoints = minPoints);
}
}
drawPrimary1 <- function(countryCode, colours = NULL, nameDisplacement = NA, resolution = 0.001,
landColour = "#FEFEE9", waterColour = "#C6ECFF", drawNames = FALSE, drawNumbers = FALSE, minPoints = 0){
if(is.na(nameDisplacement)){
if(drawNumbers){
nameDisplacement = 0;
} else {
nameDisplacement = 0.2;
}
}
require(sp);
cat("Loading Level 1 data for", countryCode, "(primary country)...", file = stderr());
con <- url(paste("http://gadm.org/data/rda/",countryCode,"_adm1.RData",sep=""));
load(con);
close(con);
cat("done!\n", file = stderr());
countryBounds <- bbox(gadm);
level1plot(gadm, colours = colours, nameDisplacement = nameDisplacement, minPoints = minPoints,
landColour = landColour, waterColour = waterColour, drawNames = drawNames, drawNumbers = drawNumbers,
resolution = resolution);
}
drawPrimary2 <- function(countryCode, colours = NULL, resolution = 0.001){
require(sp);
cat("Loading Level 2 data for", countryCode, "(primary country)...", file = stderr());
con <- url(paste("http://gadm.org/data/rda/",countryCode,"_adm2.RData",sep=""));
load(con);
close(con);
cat("done!\n", file = stderr());
level2plot(gadm, colours = colours, resolution = resolution);
}
drawPrimary3 <- function(countryCode, colours = NULL, resolution = 0.001){
require(sp);
cat("Loading Level 3 data for", countryCode, "(primary country)...", file = stderr());
con <- url(paste("http://gadm.org/data/rda/",countryCode,"_adm3.RData",sep=""));
load(con);
close(con);
cat("done!\n", file = stderr());
level3plot(gadm, colours = colours, resolution = resolution);
}
drawSecondary <- function(countryCode, landColour = "#E0E0E0", waterColour = "#C6ECFF", doLabels = FALSE,
resolution = 0.001, minPoints = 0){
require(sp);
cat("Loading Level 0 data for", countryCode, "(secondary country)...", file = stderr());
con <- url(paste("http://gadm.org/data/rda/",countryCode,"_adm0.RData",sep=""));
load(con);
close(con);
cat("done!\n", file = stderr());
countryBounds <- bbox(gadm);
level0plot(gadm, landColour = landColour, waterColour = waterColour, lineWidth = 1,
resolution = resolution, minPoints = minPoints);
if(doLabels){
plotPoint <- gadm@polygons[[1]]@labpt;
text(x = t(plotPoint), labels = as.character(gadm@data$NAME_FAO[1]));
}
}
drawCountries <- function(primaryCountry, otherCountries = NULL, drawNumbers = FALSE, doLabels = FALSE, maxLevel = 1,
resolution = 0.001, waterColour = "#C6ECFF", minPoints = 0){
require(sp);
con <- url(paste("http://gadm.org/data/rda/",primaryCountry,"_adm0.RData",sep=""));
cat("Loading Level 0 data for", primaryCountry, "(primary country)...", file = stderr());
load(con);
close(con);
primaryGADM <- gadm;
rm(gadm);
# plot image as SVG file
cat("calculating bounds... ", file = stderr());
countryBounds <- bbox(primaryGADM);
range <- apply(countryBounds,1,diff);
asp <- cos(mean(countryBounds["r2",]) * pi/180);
range[1] <- range[1] * asp;
range <- range / max(range) * 7;
cat("done!\n", file = stderr());
require(cairoDevice);
Cairo_svg(paste("GADM_",primaryCountry,"_Administrative_Divisions.svg",sep=""),
height = range[2], width = range[1]);
setupPlot(primaryGADM, waterColour = waterColour);
## draw secondary (surrounding) countries, as necessary
for(otherCode in otherCountries){
drawSecondary(otherCode, doLabels = doLabels, resolution = resolution, minPoints = minPoints);
}
if(maxLevel > 0){
drawPrimary1(primaryCountry, drawNumbers = drawNumbers,
resolution = resolution, minPoints = minPoints);
}
if(maxLevel > 1){
drawPrimary2(primaryCountry, resolution = resolution);
}
if(maxLevel > 2){
drawPrimary3(primaryCountry, resolution = resolution);
}
## draw primary country outline over the top of everything
drawPrimary0(primaryGADM, resolution = resolution, minPoints = minPoints,
landColour = NA, waterColour = NA);
dummy <- dev.off();
cat("Created '",paste("GADM_",primaryCountry,"_Administrative_Divisions.svg",sep=""),
"'\n", sep = "", file = stderr());
}
if(length(commandArgs(TRUE))>0){
## Assume first argument is the primary country, draw other countries around it
argPos <- 1;
primaryCountry <- NA;
otherCountries <- NULL;
doNumbers <- FALSE;
doLabels <- FALSE;
doLevel2 <- FALSE;
pointRes <- 0.001;
minPoints <- 0;
while(!is.na(commandArgs(TRUE)[argPos])){
argument <- commandArgs(TRUE)[argPos];
if((argument == toupper(argument)) && (nchar(argument) == 3)){
if(is.na(primaryCountry)){
primaryCountry <- argument;
cat("Adding", primaryCountry, "as primary country\n", file = stderr());
} else {
otherCountries <- c(otherCountries, argument);
cat("Adding", argument, "as secondary country\n", file = stderr());
}
} else {
if(argument == "-numbers"){
doNumbers <- TRUE;
cat("Displaying numbers at region centre points\n", file = stderr());
}
if(argument == "-minpoints"){
argPos <- argPos + 1;
minPoints <- as.numeric(commandArgs(TRUE)[argPos]);
cat("Not drawing polygons with fewer than",minPoints,"points\n", file = stderr());
}
if(argument == "-resolution"){
argPos <- argPos + 1;
pointRes <- as.numeric(commandArgs(TRUE)[argPos]);
cat("Setting point displacement resolution to", pointRes, "\n", file = stderr());
}
if(argument == "-labels"){
doLabels <- TRUE;
cat("Displaying surrounding country labels at country centre points\n", file = stderr());
}
if(argument == "-level2"){
doLevel2 <- TRUE;
cat("Will display Level 2 boundaries (if available)\n", file = stderr());
}
}
argPos <- argPos + 1;
}
maxLevel <- 1;
if(doLevel2){
maxLevel <- 2;
}
drawCountries(primaryCountry = primaryCountry, otherCountries = otherCountries, minPoints = minPoints,
drawNumbers = doNumbers, doLabels = doLabels, maxLevel = maxLevel, resolution = pointRes);
}
