Skip to content

Instantly share code, notes, and snippets.

@JoFrhwld
Created May 24, 2010 18:22
Show Gist options
  • Save JoFrhwld/412231 to your computer and use it in GitHub Desktop.
Save JoFrhwld/412231 to your computer and use it in GitHub Desktop.
Maps Functions
clip.border <- function(deldir.obj, border){
## deldir.obj = output of deldir::deldir
## border = list of borders defined in data frames
## Importantly, with columns called "x" and "y"
require(gpclib)
## bord will be a gpc polygon of all of the regions given to the argument "border"
bord <- as(matrix(ncol = 2), "gpc.poly")
for(i in seq(along = border)){
b <- as(border[[i]], "gpc.poly")
bord <- union(b, bord)
}
## Extract the polygons from the deldir object
tiles <- tile.list(deldir.obj)
out <- {}
## Keep track of polygon identity
group <- 1
for(i in seq(along = tiles)){
## Represent each tile as a gpc polygon
poly <- cbind(tiles[[i]]$x, tiles[[i]]$y)
poly <- as(poly, "gpc.poly")
## THIS IS THE CRUCIAL FUNCTION
poly <- intersect(poly, bord)
## bookkeeping and reporting
excluded <- 0
## Formatting output
if(length(poly@pts)>0){
for(j in seq(along=poly@pts)){
poly.df <- data.frame(x = poly@pts[[j]]$x, y =poly@pts[[j]]$y)
poly.df$ID <- i
poly.df$group <- group
group <- group + 1
out <- rbind(out, poly.df)
}
}else{
excluded <- excluded + 1
}
}
cat(excluded)
cat(" polygon(s) entirely outside border")
return(out)
}
## out$group defines polygon groups
## out$ID defines the data point id's
## group to ID relationship may be many to one
## out$Feature <- original.data[out$ID, ]$Feature
## Some of the data from maps is represented kind of strangely
## fix.border() is meant to reformat this data in a more useable way
## but, be sure to double check the output, since the borders may need more fixing than anticipated here.
fix.border <- function(border, rev.first = T){
out <- {}
stack <- {}
seg <- 1
for(i in 1:nrow(border)){
if(is.na(border[i,]$x)){
stack$seg <- seg
seg <- seg+1
if(is.null(out) & rev.first){
out <- rbind(out, stack[nrow(stack):1,])
}else if(is.null(out) & !rev.first){
out <- rbind(out, stack)
}else{
n.out <- nrow(out)
last <- nrow(stack)
head.dist <- sum(c((out[n.out,]$x - stack[1,]$x),(out[n.out,]$y - stack[1,]$y))^2)
tail.dist <- sum(c((out[n.out,]$x - stack[last,]$x),(out[n.out,]$y - stack[last,]$y))^2)
if(tail.dist < head.dist){
stack <- stack[nrow(stack):1,]
}
out <- rbind(out, stack[-1,])
}
stack <- {}
}else{
stack <- rbind(stack, border[i,])
}
}
out <- rbind(out, out[1,])
return(out)
}
## This is some sample code for getting boundary data from the maps package
## If the output of this data contains any NAs, you should use fix.borders(), or some other function
## usa, national
## for maps strictly of the US, or some sub region, use the "usa" map, as it is of higher quality
us.bord <- data.frame(map("usa", region = "main", plot = F)[c("x","y")])
li.bord <- data.frame(map("usa", region = "long island", plot = F)[c("x","y")])
## us.states
pa.bord <- data.frame(map("state", region = "pennsylvania", plot = F)[c("x","y")])
## some other countries
fr.bord <- data.frame(map("map", region = "France", exact = T, plot = F)[c("x","y")])
## see also, map("map", region = "France", exact = T, plot = F)$name
## This will return the borders of regions as defined by some categorical feature.
## This output is not appropriate for plotting as polygons.
## df = a data frame of tessellating polygons, as produced by clip.border()
## region = the factor defining regions
regions.border <- function(df, region){
require(gpclib)
borders <- {}
for(i in as.character(unique(df[,region]))){
## Take the subset of polygons within a region
sub <- df[df[ ,region] == i,]
poly <- as(matrix(ncol = 2), "gpc.poly")
## Take the union of all polygons
for(j in unique(sub$group)){
subsub <- subset(sub, group ==j)
subpoly <- as(subsub[,c("x","y")], "gpc.poly")
poly <- union(subpoly, poly)
}
for(j in seq(along = poly@pts)){
poly.df <- data.frame(x = poly@pts[[j]]$x, y = poly@pts[[j]]$y)
poly.df <- rbind(poly.df, poly.df[1,],c(NA,NA,NA))
borders <- rbind(borders, poly.df)
}
}
return(borders)
}
@JoFrhwld
Copy link
Author

To do:

  • Add sample plotting scripts

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment