#' gs_census
#'
#' This function takes an input landscape, a user-defined land area as an
#' [sf] polygon and returns the total population for that area for a given USA decadal census.
#' The input polygon must be projected and must overlap a USA census area
#' Data are returned for years 2000, or 2010; The default is 2010
#' For 2000 only the census block group data are available. For 2010 the block level information
#' is available but can take a long time to download for larger polygons. By default the function
#' downloads the block group level for 2010 but you can request the more detailed block data.
#' NOTE: highly recommended that you load the [tigris] package and set [options(tigris_use_cache = TRUE)]
#'
#' @param landscape A spatial polygon of class "sf" (see package [sf]). Note: the object must be projected (i.e. has a coordinate reference system) and must overlap a USA census area.
#' @param year Which census year do you want to use? In format YYYY. Currently the package supports census years 2000,
#' and 2010 (default).
#' @param spatial keep the spatial data? Default = TRUE
#' @param level Choose the census reporting level: "block_group" (default) or "block". NOTE: for year = 2000 only
#' "block_group" level is available; Both levels available for year = 2010.
#' @param api_key Your Census API key. Obtain one at http://api.census.gov/data/key_signup.html; For more information see: https://hrecht.github.io/censusapi/index.html#api-key-setup.
#'
#' @export
#' @examples
#' # NOTE: this function requires a Census API key. Obtain one at http://api.census.gov/data/key_signup.html
#' # For more information see: https://hrecht.github.io/censusapi/index.html#api-key-setup
#'
#' gs_census(buf20, year = 2010, spatial = TRUE, level = "block_group")
gs_census<-function(landscape, year = 2010, spatial = TRUE, level = "block_group",
api_key = Sys.getenv("CENSUS_API_KEY")){
# check for census api key (modified from tidycensus::get_decennial)
if (Sys.getenv("CENSUS_API_KEY") != "") {
api_key<-Sys.getenv("CENSUS_API_KEY")
} else if (is.null(api_key)) {
stop("A Census API key is required. Obtain one at http://api.census.gov/data/key_signup.html.
For more information see: https://hrecht.github.io/censusapi/index.html#api-key-setup")
}
# check that level = "block_group" for year == 2000
if (year==2000 & level == "block") stop("Block level data not available for 2000 use: level = 'block_group'")
#check that landscape is an sf polygon object
if (class(landscape)[1] != 'sf') stop("The landscape input must be an [sf] polygon object")
#check that the landscape is projected
if(is.na(sf::st_crs(landscape)$epsg) & is.na(sf::st_crs(landscape)$proj4string)) stop("The landscape object is not projected; check [st_crs]")
# check year year == 2000 or year == 2010
# load the tiger spatial data for counties
# reproject county to Albers
if(year == 2000) {
if(!exists("county2000low")) {
cty<-tigris::counties(cb=TRUE, year = 2000, class = 'sf')
cty<-sf::st_transform(cty, crs_alb)
assign("county2000low", cty, envir=globalenv())
}
county<-county2000low
} else if(year == 2010) {
if(!exists("county2010low")) {
cty<-tigris::counties(cb=TRUE, year = 2010, class = 'sf')
cty<-sf::st_transform(cty, crs_alb)
assign("county2010low", cty, envir=globalenv())
}
county<-county2010low
} else {
stop("Year must == 2000 or 2010")
}
#save input crs for output
input_crs<-sf::st_crs(landscape)
#reproject input landscape to albers (even if already in albers: just as fast)
landscape<-sf::st_transform(landscape, crs_alb) #reproject to the Albers Equal Area Conic projection used by NLCD
#get counties that landscape intersects
fips<-as.data.frame(dplyr::slice(county, unlist(sf::st_intersects(landscape, county))))
#check that fips has data. If nrow == 0 then input polygon does not overlap tiger files (i.e. not in US) or CRS wrong.
if (nrow(fips) == 0) stop("The input polygon does not overlap USA Census boundaries or projection information incorrect")
#rename fields
fips<-dplyr::select(fips, state = STATE, county = COUNTY)
#create fips string for gs_census
fips$regionin<-paste("state:", fips$state, "+county:", fips$county, sep = "")
#set up censusapi and tigris calls
if(year == 2000) {
name<-"sf1"
vars<-"P001001"
}
if(year == 2010) {
name<-"sf1"
vars<-"P0010001"
}
if(level == "block") region <- "block:*"
if(level == "block_group") region <- "block group:*"
#get the census data for total pop
pop<-list()
for(i in c(1:nrow(fips))){
pop[[i]]<-censusapi::getCensus(name = name, vintage = year, vars = vars, region = region, regionin=fips$regionin[i])
}
#rbind the pop data
tot_pop<-pop[[1]]
if(length(pop)>1) for(i in c(2:length(pop))) tot_pop<-rbind(tot_pop, pop[[i]])
names(tot_pop)[5]<-'tot_pop'
#get the census blocks/block_groups for state(s) and county(s)-each year is slightly different
#year = 2000
if(year == 2000) {
blk<-list()
for(i in c(1:nrow(fips))){
blk[[i]]<-tigris::block_groups(state = as.numeric(fips$state[i]),
county = as.numeric(fips$county[i]),year = year,
class='sf')
blk[[i]]<-dplyr::select(blk[[i]], state = STATEFP, county = COUNTYFP, tract = TRACTCE00, block.group = BLKGRPCE00,
area_land = ALAND00, area_water = AWATER00, geometry)
}
}
#year == 2010 & level == 'block_group'
if(year == 2010 & level == 'block_group') {
blk<-list()
for(i in c(1:nrow(fips))){
blk[[i]]<-tigris::block_groups(state = as.numeric(fips$state[i]),
county = as.numeric(fips$county[i]), year = year,
class='sf')
blk[[i]]<-dplyr::select(blk[[i]], state = STATEFP, county = COUNTYFP, tract = TRACTCE10, block.group = BLKGRPCE10,
area_land = ALAND10, area_water = AWATER10, geometry)
}
}
#year == 2010 & level == 'block'
if(year == 2010 & level == 'block') {
blk<-list()
for(i in c(1:nrow(fips))){
blk[[i]]<-tigris::blocks(state = as.numeric(fips$state[i]),
county = as.numeric(fips$county[i]), year = year,
class='sf')
blk[[i]]<-dplyr::select(blk[[i]], state = STATEFP, county = COUNTYFP, tract = TRACTCE10, block = BLOCKCE10,
area_land = ALAND10, area_water = AWATER10, geometry)
}
}
#rbind the blk
blocks<-blk[[1]]
if(length(blk)>1) for(i in c(2:length(blk))) blocks<-rbind(blocks,blk[[i]])
#convert area_land and area_water to numeric
blocks$area_land<-as.numeric(blocks$area_land)
blocks$area_water<-as.numeric(blocks$area_water)
#join blocks and tot_pop
blocks<-merge(blocks, tot_pop, by=names(blocks)[1:4], all.x=TRUE)
#reproject blocks
blocks<-sf::st_transform(blocks, crs_alb)
#calc area of blocks-add to attributes
blocks$blk_area<-as.numeric(sf::st_area(blocks))
#create intersection of landscape and blocks
int<-sf::st_intersection(landscape,blocks)
#calc area of intersection polygons
int$int_area<-as.numeric(sf::st_area(int))
#calc proportion of original block included in int polygons
int$pro_blk<-round(int$int_area / int$blk_area, 4)
#calc proportional population by block/block_group
int$pro_pop<-round(int$tot_pop * int$pro_blk, 2)
#prepare output list
out<-list()
#if spatial = TRUE transform intersection to input_crs and add
if(spatial) {
out$landscape<-sf::st_transform(landscape, input_crs)
out$census_raw<-sf::st_transform(blocks, input_crs)
out$census_data<-sf::st_transform(int, input_crs)
} else {
out$census_data<-as.data.frame(int)
}
#calc estimated pop for the input landscape
out$est_population<-round(sum(int$pro_pop), 2)
#calc input area in m2
out$input_area_km2<-round(as.numeric(sf::st_area(landscape) / 1000000), 2)
#calc pop density
out$est_pop_per_km2<-round(out$est_population / out$input_area_km2, 1)
#percent overlap between input landscape and census area
out$percent_overlap<-round(sum(int$int_area) / as.numeric(sf::st_area(landscape)) * 100, 1)
#add overlap warning
if(out$percent_overlap !=100) {
out$warning<-"percent_overlap not equal 100%; input landscape partially outside USA census boundaries"
warning('percent_overlap not equal 100%; input landscape partially outside USA census boundaries')
}
#return the output
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.