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