Saving Your Dorling Cartogram
You will be using your dorling cartogram file for you specific project city several times throughout the semester.
You can re-create the cartogram steps each time, or else you can simply create it once and save the new shapefile on your computer for easy access.
library( dplyr )
library( tidycensus )
library( sp )
phx.pop <-
get_acs( geography = "tract", variables = "B01003_001",
state = "AZ", county = "Maricopa", geometry = FALSE ) %>%
select( GEOID, estimate ) %>%
rename( POP=estimate )
phx.mhhi <-
get_acs( geography = "tract", variables = "B19013_001",
state = "AZ", county = "Maricopa", geometry = FALSE ) %>%
select( GEOID, estimate ) %>%
rename( MHHI=estimate )
# get a census tract shapefile
# and add census data:
library( tigris )
library( pander )
phx <- tracts( state="AZ", county="Maricopa", cb=TRUE, year=2015 )##
|
| | 0%
|
|= | 2%
|
|== | 3%
|
|==== | 5%
|
|==== | 6%
|
|===== | 8%
|
|======= | 9%
|
|======= | 10%
|
|======== | 12%
|
|========= | 12%
|
|========== | 14%
|
|=========== | 16%
|
|============= | 18%
|
|============= | 19%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 25%
|
|=================== | 27%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|======================= | 33%
|
|======================== | 34%
|
|========================= | 35%
|
|========================== | 37%
|
|=========================== | 39%
|
|============================= | 41%
|
|============================== | 43%
|
|=============================== | 44%
|
|================================ | 45%
|
|================================= | 48%
|
|=================================== | 50%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 58%
|
|========================================== | 59%
|
|========================================== | 60%
|
|============================================ | 62%
|
|============================================= | 64%
|
|============================================== | 66%
|
|================================================ | 69%
|
|================================================= | 70%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|====================================================== | 77%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================= | 87%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|======================================================================| 100%phx <- merge( phx, phx.pop, by.x="GEOID", by.y="GEOID" )
phx <- merge( phx, phx.mhhi, by.x="GEOID", by.y="GEOID" )
head( phx@data ) %>% pander()| GEOID | STATEFP | COUNTYFP | TRACTCE | AFFGEOID | |
|---|---|---|---|---|---|
| 8 | 04013040512 | 04 | 013 | 040512 | 1400000US04013040512 |
| 27 | 04013040531 | 04 | 013 | 040531 | 1400000US04013040531 |
| 42 | 04013060902 | 04 | 013 | 060902 | 1400000US04013060902 |
| 63 | 04013061027 | 04 | 013 | 061027 | 1400000US04013061027 |
| 83 | 04013061047 | 04 | 013 | 061047 | 1400000US04013061047 |
| 95 | 04013071511 | 04 | 013 | 071511 | 1400000US04013071511 |
| NAME | LSAD | ALAND | AWATER | POP | MHHI | |
|---|---|---|---|---|---|---|
| 8 | 405.12 | CT | 1988806 | 0 | 1731 | 48586 |
| 27 | 405.31 | CT | 2666330 | 6363 | 3552 | 43650 |
| 42 | 609.02 | CT | 1455458 | 0 | 3262 | 44643 |
| 63 | 610.27 | CT | 2327446 | 0 | 5454 | 66106 |
| 83 | 610.47 | CT | 2483491 | 5845 | 3863 | 119653 |
| 95 | 715.11 | CT | 2349958 | 9988 | 3674 | 58919 |
library( rgdal )
# project map and remove empty tracts
phx <- spTransform( phx, CRS("+init=epsg:3395"))
phx <- phx[ phx$POP != 0 & (! is.na( phx$POP )) , ]
# devtools::install_github( "sjewo/cartogram" )
library( cartogram ) # spatial maps w/ tract size bias reduction
library( maptools ) # spatial object manipulation
library( sf ) # 'simple features' flavor of shapefiles
# convert census tract polygons to dorling cartogram
# no idea why k=0.03 works, but it does - default is k=5
phx$pop.w <- phx$POP / 10000 # standardizes it to max of 1.5
phx_dorling <- cartogram_dorling( x=phx, weight="pop.w", k=0.03 )
plot( phx_dorling )
# install.packages( "tmap" )
library( tmap ) # thematic maps
tm_shape( phx_dorling ) +
tm_polygons( size="POP", col="MHHI", n=7, style="quantile", palette="Spectral" ) 
# WRITE TO FILE
library( geojsonio )
phx_dorling <- spTransform( phx_dorling, CRS("+proj=longlat +datum=WGS84") )
geojson_write( phx_dorling, file="phx_dorling.geojson", geometry="polygon" )## <geojson-file>
## Path: phx_dorling.geojson
## From class: SpatialPolygonsDataFrameNow you can load it from a local file, or upload it to GitHub so anyone can use the new dorling file:
library( geojsonio )
library( sp )
library( tmap )
# load from github
github.url <- "https://raw.githubusercontent.com/DS4PS/cpp-529-master/master/data/phx_dorling.geojson"
phx <- geojson_read( x=github.url, what="sp" )
plot( phx )
# from local file path
# phx <- geojson_read( "data/phx_dorling.geojson", what="sp" )
phx <- spTransform( phx, CRS("+init=epsg:3395") )
bb <- st_bbox( c( xmin = -12519146, xmax = -12421368,
ymax = 3965924, ymin = 3899074 ),
crs = st_crs("+init=epsg:3395"))
tm_shape( phx, bbox=bb ) +
tm_polygons( col="MHHI", n=10, style="quantile", palette="Spectral" ) +
tm_layout( "Dorling Cartogram", title.position=c("right","top") )
tmap_mode("view")
tm_basemap( "Stamen.Watercolor" ) +
tm_shape( phx, bbox=bb ) +
tm_polygons( col="MHHI", n=7, style="quantile", palette="-inferno" ) 