# Read data x <- read.delim("http://minato.sip21c.org/demography-special/japan-census-2010-by-pref-en.txt", stringsAsFactors=TRUE) # x <- read.delim("http://minato.sip21c.org/ldaR/japancensus2010tp.txt", encoding="CP932") # Calculating PEI (population expansion index = modified fukurami-shisuu) PEI <- function(z) { pop <- z[, 1] ag <- z[, 2] movable <- sum(pop[ag %in% c("20-24","25-29","30-34","35-39")]) fixed <- sum(pop[ag %in% c("10-14","15-19","40-44","45-49")]) return(movable/fixed*100) } # library(fmsb) # Calculate males' PEI from the data # malePEI <- as.vector(by(x$Males, x$Area, PEI, CLS=5, MODE=1)) malePEI <- as.vector(by(x[,c("Males", "Ages")], x[,"Area"], PEI)) # Classify PEI into 6 levels classes <- cut(malePEI, c(9, 9.5, 10:14)*10, include.lowest=TRUE) # or right=FALSE cols <- cm.colors(7)[-1] windows(width=1024, height=600) layout(t(1:2)) # for 2 methods of drawing map # using mapdata package (beautiful, but more difficult) if (require(maptools)==FALSE) { install.packages("maptools", dep=TRUE); library(maptools) } if (require(mapdata)==FALSE) { install.packages("mapdata", dep=TRUE); library(mapdata) } # make data xData <- data.frame( MFK=malePEI, PN=levels(x$Area) ) map("japan", type="n") for (i in xData$PN) { map("japan", region=i, fill=TRUE, add=TRUE, col=cols[classes[xData$PN==i]]) } legend("bottomright", legend=names(table(classes)), cex=1, fill=cols) title("Choropleth map for PEIs by prefecture\n (Males, Census 2010) ") # using Nippon package (easy) PN <- c("hokkaido","aomori","iwate","miyagi","akita","yamagata","fukushima", "ibaraki","tochigi","gunma","saitama","chiba","tokyo","kanagawa", "niigata","toyama","ishikawa","fukui","yamanashi","nagano","gifu", "shizuoka","aichi","mie","shiga","kyoto","osaka","hyogo","nara", "wakayama","tottori","shimane","okayama","hiroshima","yamaguchi", "tokushima","kagawa","ehime","kochi","fukuoka","saga","nagasaki", "kumamoto","oita","miyazaki","kagoshima","okinawa") classes2 <- classes[rank(PN)] if (require(NipponMap)==FALSE) { install.packages("NipponMap", dep=TRUE); library(NipponMap) } JapanPrefMap(cols[classes2], inset=FALSE) legend("bottomright", legend=names(table(classes2)), cex=1, fill=cols) title("Choropleth map for PEIs by prefecture\n (Males, Census 2010) ")