User:Gringer/gadm plot
#!/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 <[email protected]>
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);
}
Content Disclaimer
Informasi ini disarikan dari Wikipedia dan disajikan kembali untuk tujuan edukasi. Konten tersedia di bawah lisensi CC BY-SA 3.0. Kami tidak bertanggung jawab atas ketidakakuratan data yang bersumber dari kontribusi publik tersebut.
- The information displayed on this website is sourced in part or in whole from Wikipedia and has been adapted for the purpose of restating it. We strive to provide accurate and relevant information, however:
- There is no guarantee of absolute accuracy. Wikipedia is an open, collaborative project that can be edited by anyone, so information is subject to change.
- It is not intended to constitute professional advice. The content displayed is for informational and educational purposes only. For important decisions (e.g., medical, legal, or financial), please consult a professional.
- Content copyright. Wikipedia is licensed under the Creative Commons Attribution-ShareAlike License (CC BY-SA). This means that content may be reused with appropriate attribution and shared under a similar license.
- Responsible use. Any risk arising from the use of information from this website is entirely the responsibility of the user.