Preparing Village Maps of Maharashtra - Documentation
Datameet Pune
31 August 2016
Background
Datasets we have:
1. 2011 District Census handbook (DCHB) Maharashtra [xlsx and csv format]
The census has published district census handbooks available openly online which contains information on hundreds of variables (demographic, social, economic, access to water, sanitation, health, education, ICT, transport etc) for all of Maharashtra’s approximately 44000 villages.
You can go to this link and click on Maharashtra, village amenities to get the large excelsheet for all of Maharashtra (72mb) or download a csv for Maharashtra (38mb) or single district wise csv (<3mb) from this google drive link
2. Census Code Lookup Tables [csv format]
The Census codes for villages in India have not remained the same in going from census 2001 to 2011. The coding scheme was changed and lookup tables of 2001 codes and their corresponding 2011 codes are accessible on the Govt of India egov platform via this link. Note that this egov site is often problematic, alternatively you can go to this google drive link for the same file
3. Village Boundaries for Maharashtra [shapefile format and their attribute tables in csv format]
We have two versions of this shapefile and two versions of their attribute tables. Both were obtained from Bhuvan Panchayat, by accessing their WFS endpoint via using the gdal library.
Version 1: The first one has approximately 44,000 villages, it has only 2001 codes: district, sub-district, village. It does not have village names, or any 2011 codes. Hence the problem. We cannot easily append 2011 census data to this file to prepare a map.
Version 2: The second file has both 2001 and 2011 codes, district, sub-district and village and also village names. However it shows approximately 48,000 villages, many more than those existing in the 2011 DCHB.
Objectives
Study both the shapefile attribute tables
- Study version 1 and see how closely it compares with the 2011 DCHB, which can be considered the authentic source for number and names of villages
- Are there missing villages?
- Are there extra villages?
- Any missing/incorrect values for the 2001 codes
- Append 2011 codes from the census code lookup tables to the shapefile attribute table v1 and then observe
- How well do they match?
- Study version 2 and see how closely it compares with the 2011 DCHB.
- Answer questions such as are there missing villages?
- Are there extra villages?
- Any missing/incorrect values for the 2001/2011 codes?
- How is it that we find 48,000 villages instead of 44,000?
Prepare sample maps
Append the 2011 DCHB data to create thematic visualizations. These visualizations will show us spatial trends for any of the census variables we want.
Prepare a web map
More to be discussed on how this is to be done
Process
Load the required packages
library(readr)
library(data.table)
library(dplyr)
library(ggplot2)
library(VennDiagram)1. District Census Handbook (DCHB 2011)
Read in the file and View its structure
# Read it in selecting only id variables
names_dchb <-names(fread("DCHB_Village_Release_2700.csv",nrows = 0))
names_dchb_sub <-names_dchb[c(1:8,24)]
dchb <- fread("DCHB_Village_Release_2700.csv",select = c(names_dchb_sub))
# View structure of dchb
str(dchb)## Classes 'data.table' and 'data.frame': 43665 obs. of 9 variables:
## $ State Code : int 27 27 27 27 27 27 27 27 27 27 ...
## $ State Name : chr "Maharashtra" "Maharashtra" "Maharashtra" "Maharashtra" ...
## $ District Code : int 497 497 497 497 497 497 497 497 497 497 ...
## $ District Name : chr "Nandurbar" "Nandurbar" "Nandurbar" "Nandurbar" ...
## $ Sub District Code : int 3950 3950 3950 3950 3950 3950 3950 3950 3950 3950 ...
## $ Sub District Name : chr "Akkalkuwa" "Akkalkuwa" "Akkalkuwa" "Akkalkuwa" ...
## $ Village Code : int 525002 525003 525004 525005 525006 525007 525008 525009 525010 525011 ...
## $ Village Name : chr "Manibeli" "Dhankhedi" "Chimalkhadi" "Sinduri" ...
## $ Total Geographical Area (in Hectares): num 1285 768 945 1175 1181 ...
## - attr(*, ".internal.selfref")=<externalptr># dchb has 43665 observations in totalExploratory analysis: duplicates and missing values
#Any NA values in the village codes?
sum(is.na(dchb$`Village Code`))## [1] 0# Any duplicate values in village code?
sum(duplicated(dchb$`Village Code`))## [1] 0Summarizing: Count of talukas and villages in each district
# Count of total talukas and villages per district
dchb_vill_per_dist <- dchb %>% group_by(`District Name`) %>% summarise(talukas=n_distinct(`Sub District Name`),villages=n())
head(dchb_vill_per_dist)## Source: local data frame [6 x 3]
##
## District Name talukas villages
## <chr> <int> <int>
## 1 Ahmadnagar 14 1584
## 2 Akola 7 986
## 3 Amravati 14 1997
## 4 Aurangabad 9 1356
## 5 Bhandara 7 864
## 6 Bid 11 1368write.csv(dchb_vill_per_dist,"dcbh_vill_per_dist.csv",row.names = FALSE)Summarizing: Count of villages in each taluka
# Plot of total villages per district
ggplot(dchb_vill_per_dist)+geom_bar(mapping = aes(x=`District Name`, y = villages),stat = "identity") + theme(axis.text.x=element_text(angle=90))
# Count of total villages per taluka
dchb_vill_per_taluka <- dchb %>% group_by(`District Name`,`Sub District Name`) %>% summarise(Sub_district_code=unique(`Sub District Code`),villages=n())
head(dchb_vill_per_taluka)## Source: local data frame [6 x 4]
## Groups: District Name [1]
##
## District Name Sub District Name Sub_district_code villages
## <chr> <chr> <int> <int>
## 1 Ahmadnagar Akola 4201 190
## 2 Ahmadnagar Jamkhed 4214 86
## 3 Ahmadnagar Karjat 4213 120
## 4 Ahmadnagar Kopargaon 4203 80
## 5 Ahmadnagar Nagar 4209 108
## 6 Ahmadnagar Nevasa 4206 129write.csv(dchb_vill_per_taluka,"dcbh_vill_per_taluka.csv",row.names = FALSE)Summarizing: District and Sub-district, Name - Code tables (2011)
# We also create District and Sub district code tables for future reference
dstcodes<-dchb %>% group_by(`District Code`) %>% summarise(District_Name=unique(`District Name`))
head(dstcodes)## Source: local data frame [6 x 2]
##
## District Code District_Name
## <int> <chr>
## 1 497 Nandurbar
## 2 498 Dhule
## 3 499 Jalgaon
## 4 500 Buldana
## 5 501 Akola
## 6 502 Washimwrite.csv(dstcodes,"dchb_dstcode_tables.csv",row.names = FALSE)
subdstcodes<-dchb %>% group_by(`Sub District Code`) %>% summarise(Sub_District_Name = unique(`Sub District Name`))
head(subdstcodes)## Source: local data frame [6 x 2]
##
## Sub District Code Sub_District_Name
## <int> <chr>
## 1 3950 Akkalkuwa
## 2 3951 Akrani
## 3 3952 Talode
## 4 3953 Shahade
## 5 3954 Nandurbar
## 6 3955 Nawapurwrite.csv(subdstcodes,"dchb_subdstcode_tables.csv",row.names = FALSE)2. Census Code lookup (Only needed if we use Shape v1)
Read in the file and View its structure
# Read in the file
cclookup <- fread("village_2011_2001_27.csv")
# View structure of cclookup
str(cclookup)## Classes 'data.table' and 'data.frame': 43943 obs. of 10 variables:
## $ village_code_2011 : int 525002 525003 525004 525005 525006 525007 525008 525009 525010 525011 ...
## $ village_name_2011 : chr "Manibeli" "Dhankhedi" "Chimalkhadi" "Sinduri" ...
## $ state_code_2011 : int 27 27 27 27 27 27 27 27 27 27 ...
## $ district_code_2011 : int 497 497 497 497 497 497 497 497 497 497 ...
## $ sub_district_code_2011: int 3950 3950 3950 3950 3950 3950 3950 3950 3950 3950 ...
## $ village_code_2001 : int 100 200 300 400 500 600 700 800 900 1000 ...
## $ village_name_2001 : chr "Manibeli" "Dhankhedi" "Chimalkhadi" "Sinduri" ...
## $ state_code_2001 : int 27 27 27 27 27 27 27 27 27 27 ...
## $ district_code_2001 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ sub_district_code_2001: int 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, ".internal.selfref")=<externalptr># cclookup has 43943 observations in totalExploratory analysis: duplicates and missing values
#Any NA values in the village codes 2011?
sum(is.na(cclookup$village_code_2011))## [1] 0# Any duplicate values in village code 2011?
sum(duplicated(cclookup$village_code_2011))## [1] 0Exploratory analysis: villages that split up since 2001
# Find duplicate values in village code 2001 to find merged villages
cclookup_dup_2001 <- cclookup[which(duplicated(cclookup$village_code_2001)),]
# This code is incompleteSummarizing: Count of talukas and villages in each district
# Lers take a look at talukas and villages per district
cclook_vill_per_dist <- cclookup %>% group_by(district_code_2011) %>% summarise(talukas=n_distinct(sub_district_code_2011),villages=n())
head(cclook_vill_per_dist)## Source: local data frame [6 x 3]
##
## district_code_2011 talukas villages
## <int> <int> <int>
## 1 497 6 948
## 2 498 4 683
## 3 499 15 1518
## 4 500 13 1444
## 5 501 7 992
## 6 502 6 789Summarizing: Count of villages in each taluka
# Lers take a look at villages per taluka
cclook_vill_per_taluka <- cclookup %>% group_by(district_code_2011,sub_district_code_2011) %>% summarise(villages=n())
head(cclook_vill_per_taluka)## Source: local data frame [6 x 3]
## Groups: district_code_2011 [1]
##
## district_code_2011 sub_district_code_2011 villages
## <int> <int> <int>
## 1 497 3950 193
## 2 497 3951 163
## 3 497 3952 93
## 4 497 3953 184
## 5 497 3954 153
## 6 497 3955 162Comparison: CCLookup and DCHB pictorially
# Compare DCHB village codes with CCLookup village codes (2011)
# length(setdiff(dchb$`Village Code`,cclookup$village_code_2011))
# Answer: 3 villages from DCHB are missing in CCLookup
write.csv(setdiff(dchb$`Village Code`,cclookup$village_code_2011),"exclusive_dchb_missing_cclookup.csv",row.names = FALSE)
# Compare cclookup village codes with DCHB village codes (2011)
# length(setdiff(cclookup$village_code_2011,dchb$`Village Code`))
# Answer: 281 villages from CCLookup are missing in DCHB
write.csv(setdiff(cclookup$village_code_2011,dchb$`Village Code`),"exclusive_cclookup_missing_dchb.csv",row.names = FALSE)
# View the same pictorially with a Venn Diagram
grid.newpage()
draw.pairwise.venn(area1 = nrow(dchb), area2 = nrow(cclookup), cross.area = length(intersect(dchb$`Village Code`,cclookup$village_code_2011)), category = c("2011 codes DCHB","2011 Codes CCLookup"),scaled = FALSE, cat.pos = c(0,0))
## (polygon[GRID.polygon.38], polygon[GRID.polygon.39], polygon[GRID.polygon.40], polygon[GRID.polygon.41], text[GRID.text.42], text[GRID.text.43], text[GRID.text.44], text[GRID.text.45], text[GRID.text.46])Comparison: CCLookup and DCHB with a table - Taluka wise
# Does this match with DCHB 2011? We need to use the dstcode and subdstcode tables we created earlier
# Merge dstcodes lookup with cclook_vill_per_dist and save result in a comparison file
cclook_and_dchb_diff<- merge(subdstcodes,cclook_vill_per_taluka, by.x = 'Sub District Code', by.y = 'sub_district_code_2011')
cclook_and_dchb_diff <- merge(cclook_and_dchb_diff,dchb_vill_per_taluka,by.x = 'Sub District Code', by.y = 'Sub_district_code')
cclook_and_dchb_diff <- cclook_and_dchb_diff [,c(3,5,1,2,4,7)]
names(cclook_and_dchb_diff) <- c("District Code","District Name","Sub District Code","Sub District Name","Villages_cclookup","Villages_dchb")
cclook_and_dchb_diff$Diff <- cclook_and_dchb_diff$Villages_cclookup - cclook_and_dchb_diff$Villages_dchb
head (cclook_and_dchb_diff)## District Code District Name Sub District Code Sub District Name Villages_cclookup Villages_dchb
## 1 497 Nandurbar 3950 Akkalkuwa 193 190
## 2 497 Nandurbar 3951 Akrani 163 162
## 3 497 Nandurbar 3952 Talode 93 93
## 4 497 Nandurbar 3953 Shahade 184 184
## 5 497 Nandurbar 3954 Nandurbar 153 153
## 6 497 Nandurbar 3955 Nawapur 162 161
## Diff
## 1 3
## 2 1
## 3 0
## 4 0
## 5 0
## 6 1write.csv(cclook_and_dchb_diff,"cclook_and_dchb_diff.csv",row.names = FALSE)3. Shapefile Version 1
Read in the file and View its structure
# Read in the file
shapev1 <- fread("MH_Village_shapefile_attributes_v1.csv")
# View structure of shapev1
str(shapev1)## Classes 'data.table' and 'data.frame': 43897 obs. of 11 variables:
## $ gml_id : chr "village.7" "village.17" "village.24" "village.31" ...
## $ stat_name : chr "MS" "MS" "MS" "MS" ...
## $ dist_name : chr "Nashik" "Thane" "Nashik" "Nashik" ...
## $ tehs_name : chr "Yevla" "Dahanu" "Nashik" "Igatpuri" ...
## $ stcd : int 27 27 27 27 27 27 27 27 27 27 ...
## $ dstcd : int 20 21 20 20 20 20 21 21 20 21 ...
## $ sdstcd : int 15 2 11 12 12 15 2 2 10 2 ...
## $ vlgcd2011 : logi NA NA NA NA NA NA ...
## $ vlgcd2001 : int 2640700 2651400 2591600 2594200 2594000 2632400 2659400 2660500 2576600 2653300 ...
## $ shape_area: num 15515260 7289134 8875189 11413385 19658175 ...
## $ shape_len : num 16664 11242 12064 14804 21453 ...
## - attr(*, ".internal.selfref")=<externalptr># So shapev1 has 43,897 villages, but these are 2001 villages. How many villages existed in 2001 according to cclookup?Exploratory analysis: duplicates and missing values (2001)
#Any NA values in the village codes?
sum(is.na(shapev1$vlgcd2001))## [1] 579# Any duplicate values in village code?
sum(duplicated(shapev1$vlgcd2001,incomparables = NA))## [1] 205Merge Census 2011 Village Codes and Names
# Merge Census 2011 codes into shapev1 by referring to cclookup. Keep all observations from shapev1
shapev1_merged <- merge(x=shapev1,y=cclookup,by.x = 'vlgcd2001',by.y = 'village_code_2001',all.x = TRUE)
# TO BE COMPLETED!Summarizing: Count of talukas and villages in each district
Comparison: Shapev1 and DCHB pictorially
# View the overlap with DCHB pictorially with a Venn Diagram
grid.newpage()
n_12 = length(intersect(dchb$`Village Code`,cclookup$village_code_2011))
n_23 = length(intersect(cclookup$village_code_2011,unique(na.exclude(shapev1_merged$village_code_2011))))
n_13 = length(intersect(dchb$`Village Code`,unique(na.exclude(shapev1_merged$village_code_2011))))
n_123 = length(intersect(intersect(dchb$`Village Code`,cclookup$village_code_2011), unique(na.exclude(shapev1_merged$village_code_2011))))
draw.triple.venn(area1 = nrow(dchb), area2 = nrow(cclookup), area3 = length(unique(na.exclude(shapev1_merged$village_code_2011))), n12 = n_12, n23 = n_23, n13 = n_13,
n123 = n_123, category = c("DCHB", "CCLookup", "Shapev1-merged"), lty = "blank",
fill = c("blue", "green", "red"), scaled = FALSE, cat.pos = c(0,0,0), rotation.degree = 180)
## (polygon[GRID.polygon.47], polygon[GRID.polygon.48], polygon[GRID.polygon.49], polygon[GRID.polygon.50], polygon[GRID.polygon.51], polygon[GRID.polygon.52], text[GRID.text.53], text[GRID.text.54], text[GRID.text.55], text[GRID.text.56], text[GRID.text.57], text[GRID.text.58], text[GRID.text.59], text[GRID.text.60])4. Shapefile Version 2
Read in the file and View its structure
# Read in the file
shapev2 <- fread("MH_Village_shapefile_attributes_v2.csv")
shapev2$gpcode <- replace(shapev2$gpcode,which(shapev2$gpcode=="<Null>"),"")
shapev2$gpcode <- as.integer(shapev2$gpcode)
# View structure of shapev2
str(shapev2)## Classes 'data.table' and 'data.frame': 48926 obs. of 12 variables:
## $ gml_id : chr "Maharashtra_village.229" "Maharashtra_village.543" "Maharashtra_village.554" "Maharashtra_village.561" ...
## $ stname : chr "Maharashtra" "Maharashtra" "Maharashtra" "Maharashtra" ...
## $ dtname : chr "Gondia" "Gondia" "Gondia" "Gondia" ...
## $ gpname : chr "Chargaon" "Birsola" "Kasa" "Birsola" ...
## $ ipname : chr "Gondia" "Gondia" "Gondia" "Gondia" ...
## $ stcensus : int 27 27 27 27 27 27 27 27 27 27 ...
## $ dtcensus : int 507 507 507 507 507 507 507 507 507 507 ...
## $ gpcode : int 175388 175386 175418 175386 175463 175418 175415 175380 175418 175428 ...
## $ ipcode : int 4046 4046 4046 4046 4046 4046 4046 4046 4046 4046 ...
## $ villname : chr "Chargaon" "Birsola" "Kasa" "Bhadyatola" ...
## $ vlgcd2011: int 537802 537748 537747 537749 537787 537746 537750 537791 537756 537786 ...
## $ vlgcd2001: int 1276600 1271200 1271100 1271300 1275100 1271000 1271400 1275500 1272000 1275000 ...
## - attr(*, ".internal.selfref")=<externalptr>Exploratory analysis: duplicates and missing values (2011)
#Any NA values in the village codes?
sum(is.na(shapev2$vlgcd2011))## [1] 818# Any duplicate values in village code?
sum(duplicated(shapev2$vlgcd2011))## [1] 5043Summarizing: Count of talukas and villages in each district
# Do count based on unique codesCompare ShapeV2 with DCHB pictorially: Uniques only
# View the overlap with DCHB pictorially with a Venn Diagram
grid.newpage()
n_12 = length(intersect(dchb$`Village Code`,cclookup$village_code_2011))
n_23 = length(intersect(cclookup$village_code_2011,unique(na.exclude(shapev2$vlgcd2011))))
n_13 = length(intersect(dchb$`Village Code`,unique(na.exclude(shapev2$vlgcd2011))))
n_123 = length(intersect(intersect(dchb$`Village Code`,cclookup$village_code_2011), unique(na.exclude(shapev2$vlgcd2011))))
draw.triple.venn(area1 = nrow(dchb), area2 = nrow(cclookup), area3 = length(unique(na.exclude(shapev2$vlgcd2011))), n12 = n_12, n23 = n_23, n13 = n_13,
n123 = n_123, category = c("DCHB", "CCLookup", "Shapev2"), lty = "blank",
fill = c("blue", "green", "yellow"), scaled = FALSE, cat.pos = c(0,0,0), rotation.degree = 180)
## (polygon[GRID.polygon.61], polygon[GRID.polygon.62], polygon[GRID.polygon.63], polygon[GRID.polygon.64], polygon[GRID.polygon.65], polygon[GRID.polygon.66], text[GRID.text.67], text[GRID.text.68], text[GRID.text.69], text[GRID.text.70], text[GRID.text.71], text[GRID.text.72], text[GRID.text.73], text[GRID.text.74], text[GRID.text.75])Compare ShapeV2 with DCHB pictorially: Counting duplicates as well
# View the overlap with DCHB pictorially with a Venn Diagram
# Incorrect as of now
grid.newpage()
n_12 = length(intersect(dchb$`Village Code`,cclookup$village_code_2011))
n_23 = length(intersect(cclookup$village_code_2011,na.exclude(shapev2$vlgcd2011)))
n_13 = length(intersect(dchb$`Village Code`,na.exclude(shapev2$vlgcd2011)))
n_123 = length(intersect(intersect(dchb$`Village Code`,cclookup$village_code_2011), na.exclude(shapev2$vlgcd2011)))
draw.triple.venn(area1 = nrow(dchb), area2 = nrow(cclookup), area3 = length(unique(na.exclude(shapev2$vlgcd2011))), n12 = n_12, n23 = n_23, n13 = n_13,
n123 = n_123, category = c("DCHB", "CCLookup", "Shapev2"), lty = "blank",
fill = c("blue", "green", "yellow"), scaled = FALSE, cat.pos = c(0,0,0), rotation.degree = 180)
## (polygon[GRID.polygon.76], polygon[GRID.polygon.77], polygon[GRID.polygon.78], polygon[GRID.polygon.79], polygon[GRID.polygon.80], polygon[GRID.polygon.81], text[GRID.text.82], text[GRID.text.83], text[GRID.text.84], text[GRID.text.85], text[GRID.text.86], text[GRID.text.87], text[GRID.text.88], text[GRID.text.89], text[GRID.text.90])