Jump to content

User:Gringer/gadm plot

From Wikipedia, the free encyclopedia
#!/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);
}