Census Atlas Japan
library(rgdal)
library(RColorBrewer)
library("sp")
library(GISTools)
library(classInt)
library(maptools)
library("foreign")
library(gdata)
#Hacked version of the map.scale function
map.scale2 <- function (xc, yc, len, units, ndivs, subdiv = 1, tcol = "black",
scol = "black", sfcol = "black")
{
frame = par("usr")
l <- len
tic = (frame[4] - frame[3])/100
ul = l/ndivs
for (i in seq(0, ndivs - 1, by = 2)) rect(xc - l/2 + i * ul, yc, xc - l/2 + (i + 1) * ul, yc + tic/2, border = NA, col = sfcol)
lines(c(xc - l/2, xc - l/2, xc + l/2, xc + l/2), c(yc + tic, yc, yc, yc + tic), col = scol)
lines(c(xc - l/2, xc + l/2), c(yc + tic/2, yc + tic/2), col = scol)
for (i in c(0,ndivs)) text(xc - l/2 + ul * i, yc - strheight(i * subdiv) * 0.7, (i * subdiv)/100, col = tcol,cex = 0.7,font=1, family="sans")
text(xc, yc - 2 * strheight(units), units, col = tcol,cex = 0.7,font=1, family="sans")
}
#Function to deal with rounding
mround <- function(x,base){
base*round(x/base)
}
#Set working directory
setwd("/Volumes/Macintosh HD 2/Dropbox/Projects/JapanCensus/")
#setwd("/Users/alex/Dropbox/Projects/JapanCensus/")
load("JapanCensus.RData")
#Load Lookup tables
lookup <- read.xls("Cencus2010list.xls", sheet = 1)#Census lookups
muni_lookup <- read.xls("census2010municipality1901.xlsx", sheet = 1) #Municipality name lookup
muni_lookup <- muni_lookup[-1,]
muni_lookup_city <- muni_lookup[muni_lookup$Major != "",]#Subset for cities
muni_lookup <- muni_lookup[muni_lookup$Major == "",]#Remove the municipalities within cities
#Block to read the data
r_list <- list.files(path="./dbf/")
census_tables <- NULL
for (i in 1:length(r_list)){
tmp_tab_nm <- paste("CT_",gsub(".dbf","", r_list[i]),sep='')#Get table name
assign(tmp_tab_nm,read.dbf(paste("./dbf/",r_list[i],sep='')))#Read table
assign(tmp_tab_nm,get(tmp_tab_nm)[as.numeric(get(tmp_tab_nm)$areaid) > 1,c(1,12:ncol(get(tmp_tab_nm)))])#
census_tables <- c(census_tables,tmp_tab_nm)
remove(tmp_tab_nm)
}
#Get a list of the variables to be calculated
cal_vars <- lookup[!is.na(lookup$Include),c("modify","Denominator")]
#Calculate percentages for zones
for (i in 1:length(census_tables)){#table loop
#for (i in 1:length(census_tables)){#table loop
Z <- as.data.frame(get(census_tables[i])[,1])
#Next 3 commands create the area key dataframe (split up as get() and colnames() do not play nicely together! )
names(Z) <- "Geocode"; assign(paste("PCT_",census_tables[i],sep=""),Z);remove(Z) #Create a PCT table for the census table
num <- get(census_tables[i])[,2:ncol(get(census_tables[i]))]#Get census table denominators
num <- num[,intersect(colnames(num),cal_vars$modify)]#Find those numerators that we want to be converted to %
for (v in 1:ncol(num)){#loop through columns
if (all(num[v]%%1==0)){#check that variable is integer
den_id <- as.character(cal_vars[cal_vars$modify==colnames(num[v]),"Denominator"]) #Get the correct denominator
den <- get(census_tables[i])[,den_id]#Get census table numerators
tmp <- round(num[v]/den * 100,2)#Create %
assign(paste("PCT_",census_tables[i],sep=""),cbind(get(paste("PCT_",census_tables[i],sep='')),tmp))#append temp % variable onto master list
remove(tmp);remove(den_id)
}
}
remove(num);remove(den);
#print(census_tables[i])
}
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
#The following code identifies those variables with NA values - these are replaced with a zero - This replaces NaN values which are a result of 0 / 0 % calculations; also, some variables have numerator - denominator pairs that have created Inf values - i.e. the specified denominator is zero - in this instance the numerators have been assigned 0, making the assumption that the denominator is correct
#Function to look for Inf values
is.finite.data.frame <- function(obj){
sapply(obj,FUN = function(x) all(is.finite(x)))
}
#This variable is used to store the variables with Inf values
inf_var <- "test"
#Create a list of census variable tables
census_tables_PCT <- paste("PCT_",census_tables,sep="")
for (j in 1:length(census_tables_PCT)){
fix_tab <- get(census_tables_PCT[j])#Create a temp copy of census table
fix_tab[is.na(fix_tab)] <- as.numeric(0)#change NA to zero
inf_var_tmp <- colnames(fix_tab)[!is.finite.data.frame(fix_tab)]#This gets a list of columns with Inf values
inf_var <- c(inf_var,inf_var_tmp)#Add the colums with Inf values to the master list
is.na(fix_tab) <- do.call(cbind,lapply(fix_tab, is.infinite))#Change Inf values to NA
fix_tab[is.na(fix_tab)] <- as.numeric(0)#change NA to zero
assign(census_tables_PCT[j],fix_tab)#Overwrite the PCT table
}
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
#save.image("JapanCensus.RData")
# # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # #
#Create Plots # # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # #
# # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # #
#setwd("/Users/alex/Dropbox/Projects/JapanCensus/")
setwd("/Volumes/Macintosh HD 2/Dropbox/Projects/JapanCensus/")
#Read Shapefile Zones
#zone <- readOGR(".", "census2010_Zone")
#zone <- zone[-1,]
zone <- readOGR(".", "census2010_chocho_projected_nosea_nolake_final")
#Read Shapefile Places, then append a population count for the zone
places <- readOGR("gazateer_japan", "JP_Place_Names_Muni_Projected")
places <- places[!is.na(places@data$AREA_CODE),]
zone_places <- zone
zone_places@data <- data.frame(zone_places@data, CT_2010_002a[match(zone_places@data[, "hgeocode"], CT_2010_002a[, "Geocode"]), ])#append data to map for municipality
zone_places@data <- zone_places@data[,c("hgeocode","to6_2_0000")]
o <- over(places, zone_places)
o <- o[,"to6_2_0000"]
places@data <- cbind(places@data,o)
remove(zone_places);remove(o)
#Create a list containing the municipality codes
#municipality <- as.character(unique(substr(zone@data$hgeocode,1,5)))
municipality <- as.character(muni_lookup$AREA_CODE)
#setwd("/Users/alex/Dropbox/Projects/JapanCensus/maps/")
#setwd("/Volumes/Macintosh HD 2/Dropbox/Projects/JapanCensus/maps/")
setwd("/Users/alex/Desktop/map_junk")
#Create plot
for (i in 9:length(municipality)){#Munipality loop
print(c(municipality[i],i))
for (j in 1:length(census_tables_PCT)){#Census table loop
#Read and join data
tmp_zone <- zone[substr(zone@data$hgeocode,1,5) == municipality[i],] #subset of the polygons for the municipality
tmp_zone@data = data.frame(tmp_zone@data, get(census_tables_PCT[j])[match(tmp_zone@data[, "hgeocode"], get(census_tables_PCT[j])[, "Geocode"]), ])#append data to map for municipality
tmp_place <- places[places@data$AREA_CODE == municipality[i],]#get places within municipality
tmp_place <- tmp_place[order(tmp_place$o,decreasing = TRUE),]#order by population
#loop to check the number of areas identified in geobase
if (nrow(tmp_place) >= 5) {
if (nrow(tmp_place) > 14) {
tmp_place <- tmp_place[1:15,] #Select random 5 points from top 15
tmp_place <- tmp_place[sample(1:15, 5, replace=F),]
} else {
tmp_place <- tmp_place[sample(1:nrow(tmp_place), 5, replace=F),]#Select random 5 points
}
}
#get variables to plot
vars_to_plot <- colnames(tmp_zone@data[(grep("Geocode", colnames(tmp_zone@data))+1):ncol(tmp_zone@data)])
for (k in 1:length(vars_to_plot)){#Start variable loop
#Get breaks
tmp_plot_vars <- tmp_zone@data[,paste(vars_to_plot[k])]
tmp_plot_vars[is.nan(tmp_plot_vars)] <- as.numeric(0)#This replaces NaN values which are a result of 0 / 0 % calculations
tmp_plot_vars <- tmp_plot_vars[!is.nan(tmp_plot_vars)]
if (length(unique(tmp_plot_vars)) > 5) {#Loop to check if there are more than 5 unique values
#Create colour pallet
#Loop to check if there are less than 5 zones
if (nrow(tmp_zone) >=5) {
my_colours <- brewer.pal(5, "Blues")
breaks <- classIntervals(tmp_plot_vars, n = 5, style = "fisher", unique = TRUE)
breaks <- breaks$brks
} else {
my_colours <- brewer.pal((nrow(tmp_zone)), "Blues")
breaks <- classIntervals(tmp_plot_vars, n = (nrow(tmp_zone)), style = "fisher", unique = TRUE)
breaks <- breaks$brks
}
#Create map output as PDF
pdf(paste("MAP_",municipality[i],"_VAR___",census_tables_PCT[j],"___",vars_to_plot[k],".pdf",sep=''))
par(mfrow=c(2,1), oma=c(0,0,0,0),mar=c(0, 0, 0, 0),bty="n",family="serif",bg = '#FFFFFF')
plot(tmp_zone, col = my_colours[findInterval(tmp_zone@data[,vars_to_plot[k]], breaks, all.inside = TRUE)], axes = FALSE, border = NA)
plot(tmp_place,add=TRUE,pch=19,cex=.5,col="red")#Plots 5 places
pointLabel(tmp_place@coords[,"coords.x1"],tmp_place@coords[,"coords.x2"], labels = tmp_place@data[,"Name"], cex = 0.7,font=1, family="sans")
plot(tmp_zone, col = NA, axes = FALSE, border = NA)
#work out scale position
plot_width <- (tmp_zone@bbox[3] - tmp_zone@bbox[1])/1000
scale_width <- mround(plot_width,10) *100
if (scale_width == 0) { scale_width <- 1000}
scale_break <- 500
scale_divs <- scale_width / scale_break
X <- tmp_zone@bbox[3] - (scale_width * 1.2)
Y <- tmp_zone@bbox[4] - (scale_width * 1.1)
map.scale2(X,Y,scale_width,"km",scale_divs,scale_break,sfcol='#08519C')
dev.off()
#Create legend output as PDF
pdf(paste("LEG_",municipality[i],"_VAR___",census_tables_PCT[j],"___",vars_to_plot[k],".pdf",sep=''), family="sans")
plot.new()
legend("topleft", legend = leglabs(round(breaks, digits = 1), between = " to "), fill = my_colours, bty = "o", ,bg="#FFFFFF", border = NA, title="Percent" ,box.lty=0)
dev.off()
#Trim the PDF that were created to remove whitespace around the plots - this uses an external library that runs on OSX (need PDFCrop installed - http://users.skynet.be/tools/PDFTools.tgz installed)
#crop map
system(paste("pdfcrop '",getwd(),"/",paste("MAP_",municipality[i],"_VAR___",census_tables_PCT[j],"___",vars_to_plot[k],".pdf'",sep='')," '",getwd(),"/",paste("MAP_",municipality[i],"_VAR___",census_tables_PCT[j],"___",vars_to_plot[k],".pdf'",sep=''),sep=''),wait=TRUE,ignore.stdout = TRUE, ignore.stderr = TRUE)
#crop legend
system(paste("pdfcrop '",getwd(),"/",paste("LEG_",municipality[i],"_VAR___",census_tables_PCT[j],"___",vars_to_plot[k],".pdf'",sep='')," '",getwd(),"/",paste("LEG_",municipality[i],"_VAR___",census_tables_PCT[j],"___",vars_to_plot[k],".pdf'",sep=''),sep=''),wait=TRUE,ignore.stdout = TRUE, ignore.stderr = TRUE)
}#Loop to check if the values are all 0
else{
next
}
}#Variable loop
summary(table(tmp_zone@data[,vars_to_plot[k]]))
remove(tmp_zone)
}#Table loop
# # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # #
#Make Maps # # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # #
# # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # ## # #
#This chunk creates the latex code for displaying the maps
#setwd("/Users/alex/Dropbox/Projects/JapanCensus/maps/")
#setwd("/Volumes/Macintosh HD 2/Dropbox/Projects/JapanCensus/maps/")
#Create a list of the maps, and a table - variable lookup
maps_list <- list.files(path=".")
if (length(maps_list) > 0) {#check that there are some images generated
maps_list <- maps_list[1:length(maps_list)-1]#remove .Rnw from the list
maps_list <- gsub("LEG_","",maps_list)
maps_list <- gsub("MAP_","",maps_list)
maps_list <- gsub(".pdf","",maps_list)
maps_list <- gsub("_VAR","",maps_list)
maps_list <- gsub("PCT_","",maps_list)
maps_list <- strsplit(maps_list,"___")
maps_list <- data.frame(matrix(unlist(maps_list), nrow=length(maps_list), byrow=T))
colnames(maps_list) <- c("Muni","C_Tab","C_Var")
maps_list <- unique(maps_list)
#muni_list <- as.character(unique(maps_list$Muni))#Create unique municipalities list
muni_list <- municipality[i]
for (i in 1:length(muni_list)){#loop for each municipality
#Create a municipality tex file
file.create(paste(muni_list[i],".tex",sep=''))
#Open tex file for edits
fileConn<-file(paste(muni_list[i],".tex",sep=''))
#Get the names of the municipality
muni_name <- as.character(muni_lookup[muni_lookup$AREA_CODE == muni_list[i],"NAME_ENG"])
header <- paste("\\documentclass[a4paper,11pt]{article}
\\pagenumbering{gobble}
\\makeatletter
\\renewcommand{\\@dotsep}{10000}
\\makeatother
\\usepackage{helvet}
\\renewcommand{\\familydefault}{\\sfdefault}
\\setlength{\\textwidth}{15cm}
\\setlength{\\oddsidemargin}{1cm}
\\setlength{\\evensidemargin}{1cm}
\\setlength{\\topmargin}{1cm}
\\title{",paste("Open Census Atlas: Japan 2010 (",muni_name," / ", muni_list[i],")",sep=''),"}
\\author{Alex D Singleton}
\\usepackage[hidelinks]{hyperref}
\\usepackage{float}
\\usepackage[pdftex]{graphicx}
\\usepackage[margin=2cm]{geometry}
\\usepackage{subfigure}
\\usepackage{caption}
\\usepackage{graphicx}
\\begin{document}
\\maketitle
\\label{listfigs}
\\listoffigures
\\clearpage",sep="\n")
#Subset a list of tables / variables for selected municipality
muni_temp <- maps_list[maps_list$Muni == muni_list[i],]
muni_temp_tab_list <- as.character(unique(muni_temp$C_Tab))
#Some initial variables for the figures loop
##content_counter <- 1
fig_list <- NULL
for (t in 1:length(muni_temp_tab_list)){#loop through tables within municipality
muni_temp_tab_var_list <- as.character(unique(muni_temp[muni_temp$C_Tab == muni_temp_tab_list[t],"C_Var"]))#list of variables
for (v in 1:length(muni_temp_tab_var_list)){#loop through variables within tables within municipality
table_name <- as.character(unique(lookup[lookup$modify == muni_temp_tab_var_list[v],"Census2010.contents"]))#get table name
map_id <- paste(muni_list[i],"_VAR___PCT_",muni_temp_tab_list[t],"___",muni_temp_tab_var_list[v],sep="")#create the map id
variable_name <- as.character(unique(lookup[lookup$modify == muni_temp_tab_var_list[v],"Numerator_Label"]))#get variable name
##assign(paste("content_",content_counter,sep=''), paste("\\begin{figure}","\\hyperref[listfigs]{\\includegraphics[width=\\textwidth]{MAP_",map_id,".pdf}}\n","\\includegraphics[width=3cm]{LEG_",map_id,".pdf}\n","\\label{fig:",muni_temp_tab_var_list[v],"}\n","\\caption{",paste(table_name," [",variable_name,"]",sep=''),"}\\end{figure} \\clearpage",sep = ""))
#assign(paste("content_",muni_temp_tab_var_list[v],sep=''), paste("\\begin{figure}\\hyperref[listfigs]{\\includegraphics[width=15cm]{MAP_",map_id,".pdf}}\n\\includegraphics[width=3cm]{LEG_",map_id,".pdf}\n","\\label{fig:",muni_temp_tab_var_list[v],"}\n","\\caption{",paste(table_name," [",variable_name,"]",sep=''),"}\\end{figure} \\clearpage",sep = ""))
#assign(paste("content_",muni_temp_tab_var_list[v],sep=''), paste("\\begin{center}\\begin{figure}\\hyperref[listfigs]{\\includegraphics[width=15cm,height=20cm,keepaspectratio]{MAP_",map_id,".pdf}}\\end{center}\\\\\\includegraphics[width=3cm]{LEG_",map_id,".pdf}\n","\\label{fig:",muni_temp_tab_var_list[v],"}\n","\\caption{",paste(table_name," [",variable_name,"]",sep=''),"}\\end{figure} \\clearpage",sep = ""))
assign(paste("content_",muni_temp_tab_var_list[v],sep=''),paste("
\\begin{figure}[H]
\\centering
\\hyperref[listfigs]{\\includegraphics[width=15cm,height=20cm,keepaspectratio]{MAP_",map_id,".pdf}}
\\end{figure}
\\begin{figure}[H]
\\includegraphics[width=3cm]{LEG_",map_id,".pdf}
\\caption{", paste(table_name," [",variable_name,"]",sep=''),"}
\\end{figure}\\clearpage",sep=""))
##fig_list <- c(fig_list,get(paste("content_",content_counter,sep='')))
fig_list <- c(fig_list,get(paste("content_",muni_temp_tab_var_list[v],sep='')))
#remove(list=paste("content_",content_counter,sep=''))
remove(list=paste("content_",muni_temp_tab_var_list[v],sep=''))
remove(map_id)
remove(variable_name)
remove(table_name)
##content_counter <- content_counter + 1
}#End variable loop
}#End table loop
#Adds content to the tex file
footer <- "\\end{document}"
content <- as.character(fig_list)
writeLines(c(header,content,footer), fileConn)
close(fileConn)
#Create LaTex document
system(paste("pdflatex '",getwd(),"/",muni_list[i],".tex'",sep=''),wait=TRUE,ignore.stdout = TRUE, ignore.stderr = TRUE)
system(paste("pdflatex '",getwd(),"/",muni_list[i],".tex'",sep=''),wait=TRUE,ignore.stdout = TRUE, ignore.stderr = TRUE)
del_ext <- c(".aux",".lof",".log",".out")
file.remove(paste(muni_list[i],del_ext,sep=''))
#Adds a cover (Uses pdftk - package versions downloaded from http://www.pdflabs.com/tools/pdftk-the-pdf-toolkit/)
#cover_location <- "/Users/alex/Dropbox/Projects/JapanCensus/Cover.pdf"
cover_location <- "'/Volumes/Macintosh HD 2/Dropbox/Projects/JapanCensus/Cover.pdf'"
atlas_location <- paste("'",getwd(),"/",muni_list[i],".pdf'",sep='')
#final_atlas_location <- paste("'/Volumes/cloudshare/share/Alex/",muni_list[i],".pdf'",sep='')
final_atlas_location <- paste("'/Volumes/ALEX/ATLAS/",muni_list[i],".pdf'",sep='')
#final_atlas_location <- paste("'/Volumes/Macintosh HD 2/Dropbox/Projects/JapanCensus/ATLAS/",muni_list[i],".pdf'",sep='')
system(paste("pdftk",cover_location,atlas_location,"cat output",final_atlas_location),wait = TRUE)
file.remove(paste(muni_list[i],".pdf",sep=''))
system(paste("find ",paste("'",getwd(),"'",sep='')," -name '*.pdf' -delete",sep=''),wait = TRUE)
system(paste("find ",paste("'",getwd(),"'",sep='')," -name '*.tex' -delete",sep=''),wait = TRUE)
}#End municipality loop - latex
}#check if images exist
}#Municipality loop - map generation