24 Days of R: Day 21

I drove through Tennessee today. Despite having grown up in Kentucky, I had never really noticed the Tennessee flag, but earlier this year, while in Nashville, I saw it and thought it was one of the cooler state flags. I downloaded a PNG file of the flag from wikicommons and tried to do what I thought would be a very simple thing. I wanted to draw the state of Tennessee and have the flag appear on top, but clipped along the borders of the state.

Drawing a picture of the flag is fairly easy.

library(raster)
library(png)
pngTN = readPNG("./figure/800px-Tennessee_state_flag.png")
plot(x = c(0, 1000), y = c(0, 1000), type = "n", xaxt = "n", yaxt = "n", xlab = "", 
    ylab = "", bty = "n")
rasterImage(pngTN, 0, 0, 1000, 1000)

plot of chunk TN_Flag

As is drawing a picture of the state.

library(maps)
tn = getData("GADM", country = "USA", level = 1)
tn = tn[match("tennessee", tolower(tn$NAME_1)), ]
plot(tn)

plot of chunk TN

I can even draw one on top of the other.

tn.bbox = tn@bbox
xleft = tn.bbox[1, 1]
xright = tn.bbox[1, 2]
ybottom = tn.bbox[2, 1]
ytop = tn.bbox[2, 2]
plot(x = c(xleft, xright), y = c(ybottom, ytop), type = "n", xaxt = "n", yaxt = "n", 
    xlab = "", ylab = "", bty = "n")
rasterImage(pngTN, xleft, ybottom, xright, ytop)
plot(tn, add = TRUE)

plot of chunk Both

This next bit required about an hour of work, fussing and research and it still doesn't work right. The assignment of the raster values to the RasterLayer won't accept character data. I converted to factors, which means that the shape of the image is fine, but the colors are off. It's also a very slow operation.

library(raster)
library(rgdal)
mojo = raster(xmn = xleft, xmx = xright, ymn = ybottom, ymx = ytop, nrow = 480, 
    ncol = 800)
rasterTN = as.raster(pngTN)
rasterTN = as.factor(rasterTN)
mojo[] = rasterTN
masked = mask(mojo, tn)
## Found 1 region(s) and 1 polygon(s)
plot(masked)

plot of chunk Clipped

Until tonight, I had thought of a “raster” as something that sounds like a person from Jamaica. It's clear that I need to learn more as this is a hot topic in GIS circles.

As usual, there's virtually nothing new here. Here are a few links to other dabblings, inspiration and full code examples.:
Use of mask is something that I saw on gis.stackexchange

I came across this site while doing basic interweb research. I didn't spend a great deal of time with that code as it was a bit mysterious.

Last year, is.R() did a number of posts about flags. This is where I first saw the use of the rasterImage function

sessionInfo
## function (package = NULL) 
## {
##     z <- list()
##     z$R.version <- R.Version()
##     z$platform <- z$R.version$platform
##     if (nzchar(.Platform$r_arch)) 
##         z$platform <- paste(z$platform, .Platform$r_arch, sep = "/")
##     z$platform <- paste0(z$platform, " (", 8 * .Machine$sizeof.pointer, 
##         "-bit)")
##     z$locale <- Sys.getlocale()
##     if (is.null(package)) {
##         package <- grep("^package:", search(), value = TRUE)
##         keep <- sapply(package, function(x) x == "package:base" || 
##             !is.null(attr(as.environment(x), "path")))
##         package <- sub("^package:", "", package[keep])
##     }
##     pkgDesc <- lapply(package, packageDescription, encoding = NA)
##     if (length(package) == 0) 
##         stop("no valid packages were specified")
##     basePkgs <- sapply(pkgDesc, function(x) !is.null(x$Priority) && 
##         x$Priority == "base")
##     z$basePkgs <- package[basePkgs]
##     if (any(!basePkgs)) {
##         z$otherPkgs <- pkgDesc[!basePkgs]
##         names(z$otherPkgs) <- package[!basePkgs]
##     }
##     loadedOnly <- loadedNamespaces()
##     loadedOnly <- loadedOnly[!(loadedOnly %in% package)]
##     if (length(loadedOnly)) {
##         names(loadedOnly) <- loadedOnly
##         pkgDesc <- c(pkgDesc, lapply(loadedOnly, packageDescription))
##         z$loadedOnly <- pkgDesc[loadedOnly]
##     }
##     class(z) <- "sessionInfo"
##     z
## }
## <bytecode: 0x0000000013b2c778>
## <environment: namespace:utils>
Advertisements
Posted in R

2 thoughts on “24 Days of R: Day 21

  1. Easier way is to read the PNG into a raster brick rather than a PNG, then assign it an extent, then use mask….

    flag = brick(“./800px-Tennessee_state_flag.png”)
    extent(flag)=extent(tn)
    plotRGB(mask(flag,tn))
    plot(tn,add=TRUE,lwd=2) # add lines

    Now the problem here is that the flag is stretched because the aspect ratio is wrong. A bit of trial and error gives me this:

    extent(flag)=extent(-90.3,xmax=-81.6,ymin=33.5,ymax=38)
    plotRGB(mask(flag,tn))
    plot(tn,add=TRUE,lwd=2) # add lines

    which keeps the circle circular.

    Only 49 more states to go!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s