Maps Functions

clip.border.R

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

fix.border.R

## 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)	
}

getting.maps.R

## 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

regions.border.R

## 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)
}

To do:

• Add sample plotting scripts