Most of us end up having to work with Census data on a regular basis. In the past, we used American Factfinder (replaced by the data.census.gov) or then the Interuniversity Consortium of Political and Social Research (ICPSR) databases to grab the data we needed. Mercifully, with the opening up of a lot of government data, accessing Census data has become very easy. There are several packages that allow you to do so but I will initially focus on two packages – {tidycensus} and {censusapi} – given their ease of use. In addition to these packages, I will show you how to work with APIs to get data from the World Bank, the U.S. Geological Survey (USGS), and other sources. Please also be sure to look at IPUMS
Most of us end up having to work with Census data on a regular basis. In the past, we used American Factfinder or then the ICPSR databases to grab the data we needed. Mercifully, with the opening up of a lot of government data grabbing Census data has become very easy. There are several packages that allow you to do so but I will focus on two packages – {tidycensus} and {censusapi} – given their ease of use. In addition to these packages, I will also show you how to use APIs to access data from the World Bank and from the U.S. Geological Survey (USGS).
{tidycensus}Let us start with {tidycensus}. As usual, we install the package and, since we will need it, get a Census API key. You can signup for the key here. Check your email and write your key somewhere and remember it. Install {tidycensus} next.
library(tidycensus)
library(tidyverse)
census_api_key("YOUR API KEY GOES HERE", install = TRUE)Your key will be saved to .Renviron and automatically used in future sessions with {tidycensus} so be sure to not use the install = TRUE switch for subsequent runs otherwise you will be asked every time if you want to overwrite the previous setting in .Renviron.
The two data-streams of interest will either be the decennial census or then, most likely, the regularly updated American Community Survey (ACS) series. You can get each with specific commands:
load_variables(
year = 2016,
dataset = "acs5",
cache = TRUE
) -> my.vars.acs
load_variables(
year = 2000,
dataset = "sf3",
cache = TRUE
) -> my.vars.d.sf3 What if I am interested in the Summary File 1 from 2010 or Summary File 3 from 2000?
load_variables(
year = 2010,
dataset = "sf1",
cache = TRUE
) -> my.vars.d.sf1
load_variables(
year = 2000,
dataset = "sf3",
cache = TRUE
) -> my.vars.d.sf3 Each of these will show you what variables are available, and make it easier for you to choose the ones you want to work with. Assume we want total population, variable P0010001, from the decennial census
get_decennial(
geography = "state",
variables = c(totpopn = "P001001"),
year = 2000
) -> state.popn.d
get_decennial(
geography = "county",
variables = c(totpopn = "P001001"),
state = "OH",
year = 2000
) -> county.popn.d
get_decennial(
geography = "tract",
variables = c(totpopn = "P001001"),
state = "OH",
county = "Athens",
year = 2000
) -> tract.popn.d and then from the ACS. Note that the default for each is 2020 and the most recent ACS (2019-2023).
get_acs(
geography = "state",
variables = c(totpopn = "B01003_001"),
year = 2023
) -> state.popn.acs
get_acs(
geography = "county",
variables = c(totpopn = "B01003_001"),
state = "OH",
year = 2023
) -> county.popn.acs
get_acs(
geography = "tract",
variables = c(totpopn = "B01003_001"),
state = "OH",
county = "Athens",
year = 2023
) -> tract.popn.acs You can also get an entire table instead of having to list variables, such as shown below, with the default decennial census for the package’s current version (2010).
get_decennial(
geography = "county",
table = c("P012"),
year = 2000
) -> county.table.d
get_decennial(
geography = "state",
table = c("P012"),
year = 2000
) -> state.table.d Population data for all tracts in the country?
library(tidycensus)
library(purrr)
unique(fips_codes$state)[1:51] -> us
map_df(us, function(x) {
get_acs(
geography = "tract",
variables = "B01003_001",
state = x,
year = 2023,
geometry = TRUE
)
}
) -> totalpop Often you want to convert some measure into a percent. You can do that either by also downloading the denominator, or then leaning on built-in functions (as shown below). In this example, I am grabbing the B02001 table for all tracts in Cuyahoga county, Ohio.
get_acs(
"tract",
county = "Cuyahoga",
state = "Ohio",
year = 2023,
table = "B02001",
summary_var = "B02001_001",
geometry = TRUE
) -> cuy | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 5% | |=== | 6% | |==== | 7% | |===== | 8% | |===== | 9% | |====== | 10% | |======= | 11% | |======== | 13% | |======== | 14% | |========= | 15% | |========== | 16% | |========== | 17% | |=========== | 19% | |============ | 20% | |============= | 21% | |============== | 23% | |============== | 24% | |=============== | 25% | |================ | 26% | |================ | 27% | |================= | 28% | |================== | 29% | |================== | 30% | |=================== | 32% | |==================== | 33% | |==================== | 34% | |===================== | 35% | |====================== | 36% | |====================== | 37% | |======================= | 38% | |======================== | 39% | |======================== | 41% | |========================= | 42% | |========================== | 43% | |========================== | 44% | |=========================== | 45% | |============================ | 46% | |============================ | 47% | |============================= | 49% | |============================== | 50% | |=============================== | 51% | |=============================== | 52% | |================================ | 53% | |================================= | 55% | |================================== | 56% | |=================================== | 58% | |=================================== | 59% | |==================================== | 60% | |===================================== | 61% | |===================================== | 62% | |====================================== | 63% | |======================================= | 64% | |======================================= | 66% | |======================================== | 67% | |========================================= | 68% | |========================================= | 69% | |========================================== | 70% | |=========================================== | 71% | |============================================ | 73% | |============================================= | 75% | |============================================= | 76% | |============================================== | 77% | |=============================================== | 78% | |=============================================== | 79% | |================================================ | 80% | |================================================= | 81% | |================================================== | 83% | |================================================== | 84% | |=================================================== | 85% | |==================================================== | 86% | |==================================================== | 87% | |===================================================== | 88% | |====================================================== | 89% | |======================================================= | 92% | |======================================================== | 93% | |======================================================== | 94% | |========================================================= | 95% | |========================================================== | 96% | |========================================================== | 97% | |=========================================================== | 98% | |============================================================| 99% | |============================================================| 100%cuy %>%
head()Simple feature collection with 6 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.71021 ymin: 41.41891 xmax: -81.69512 ymax: 41.4326
Geodetic CRS: NAD83
GEOID NAME variable
1 39035105900 Census Tract 1059; Cuyahoga County; Ohio B02001_001
2 39035105900 Census Tract 1059; Cuyahoga County; Ohio B02001_002
3 39035105900 Census Tract 1059; Cuyahoga County; Ohio B02001_003
4 39035105900 Census Tract 1059; Cuyahoga County; Ohio B02001_004
5 39035105900 Census Tract 1059; Cuyahoga County; Ohio B02001_005
6 39035105900 Census Tract 1059; Cuyahoga County; Ohio B02001_006
estimate moe summary_est summary_moe geometry
1 2584 447 2584 447 MULTIPOLYGON (((-81.70976 4...
2 2126 411 2584 447 MULTIPOLYGON (((-81.70976 4...
3 138 90 2584 447 MULTIPOLYGON (((-81.70976 4...
4 0 13 2584 447 MULTIPOLYGON (((-81.70976 4...
5 0 13 2584 447 MULTIPOLYGON (((-81.70976 4...
6 0 13 2584 447 MULTIPOLYGON (((-81.70976 4...Now I can run mutate(...) to convert each estimate into its corresponding percentage.
cuy %>%
mutate(
pct = round(estimate / summary_est) * 100
) -> cuy
cuy %>%
head()Simple feature collection with 6 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.60645 ymin: 41.52275 xmax: -81.5981 ymax: 41.53004
Geodetic CRS: NAD83
GEOID NAME variable
1 39035118400 Census Tract 1184, Cuyahoga County, Ohio B02001_001
2 39035118400 Census Tract 1184, Cuyahoga County, Ohio B02001_002
3 39035118400 Census Tract 1184, Cuyahoga County, Ohio B02001_003
4 39035118400 Census Tract 1184, Cuyahoga County, Ohio B02001_004
5 39035118400 Census Tract 1184, Cuyahoga County, Ohio B02001_005
6 39035118400 Census Tract 1184, Cuyahoga County, Ohio B02001_006
estimate moe summary_est summary_moe geometry
1 1203 199 1203 199 MULTIPOLYGON (((-81.60645 4...
2 18 18 1203 199 MULTIPOLYGON (((-81.60645 4...
3 1137 193 1203 199 MULTIPOLYGON (((-81.60645 4...
4 11 17 1203 199 MULTIPOLYGON (((-81.60645 4...
5 0 11 1203 199 MULTIPOLYGON (((-81.60645 4...
6 0 11 1203 199 MULTIPOLYGON (((-81.60645 4...
pct
1 100
2 0
3 100
4 0
5 0
6 0The White count is B02001_002 and Black is B02001_003
Say I wanted to map just these two for the tracts.
cuy %>%
filter(
variable %in% c('B02001_002', 'B02001_003')
) %>%
mutate(
popgroup = case_when(
variable == "B02001_002" ~ 'White Alone',
variable == "B02001_003" ~ 'Black Alone'
)
) -> cuy_df
cuy_df %>%
ggplot() +
geom_sf(
aes(fill = pct)
) +
facet_wrap(~ popgroup) +
theme_void() +
theme(legend.position = 'top') +
scale_fill_viridis_c(option = "rocket") +
labs(
fill = "Percent"
)I threw in the mapping here to show how easy it would be to pipe through from getting data to visualizing it.
{censusapi}This package will allow you to grab a vast array of Census products, and very easily I might add.
library(censusapi)
listCensusApis() -> apis_df Application Programming Interfaces (APIs) could have slightly varying parameters so it pays to check the API documentation available here.
It is easy to find variable names via the built-in function listCensusMetadata. As an example, the bureau’s small area health insurance estimates are shown below, as well as the small area income and poverty estimates.
listCensusMetadata(
name = "timeseries/healthins/sahie",
type = "variables"
) -> sahie_vars
listCensusMetadata(
name = "timeseries/poverty/saipe",
type = "variables"
) -> saipe_vars Curious about available geographies? Switch type = to geography:
listCensusMetadata(
name = "timeseries/healthins/sahie",
type = "geography"
) -> sahie_geos
listCensusMetadata(
name = "timeseries/poverty/saipe",
type = "geography"
) -> saipe_geos Grab the most recent county-level data but note that the latest SAHIE are for 2015 while the latest SAIPE are for 2016. Every variable is downloaded as a chr so you will need to flip what should be numeric into numeric.
getCensus(
name = "timeseries/healthins/sahie",
vars = c("NAME", "IPRCAT", "IPR_DESC", "PCTUI_PT"),
region = "county:*",
regionin = "state:39",
time = 2019
) -> sahie_counties
head(sahie_counties, n = 12L)
getCensus(
name = "timeseries/poverty/saipe",
vars = c("NAME", "SAEPOVRTALL_PT", "SAEMHI_PT"),
region = "county:*",
regionin = "state:39",
time = 2019
) -> saipe_counties
head(saipe_counties, n = 12L)
saipe_counties %>%
mutate(
prate = as.numeric(SAEPOVRTALL_PT),
mdhinc = as.numeric(SAEMHI_PT)
) -> saipe_countiesIf you want data for a lot of geographies, the package will let you do it in seconds. For example, if you want tract-level data for all tracts, you can get it as shown below:
fips [1] "01" "02" "04" "05" "06" "08" "09" "10" "11" "12" "13" "15" "16"
[14] "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29"
[27] "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42"
[40] "44" "45" "46" "47" "48" "49" "50" "51" "53" "54" "55" "56"tracts <- NULL
for (f in fips) {
stateget <- paste("state:", f, sep = "")
temp <- getCensus(
name = "dec/sf1",
vintage = 2010,
vars = c("NAME", "P001001", "H010001"),
region = "tract:*",
regionin = stateget
)
tracts <- rbind(tracts, temp)
}
head(tracts) state county tract NAME
1 01 001 020100 Census Tract 201, Autauga County, Alabama
2 01 001 020500 Census Tract 205, Autauga County, Alabama
3 01 001 020300 Census Tract 203, Autauga County, Alabama
4 01 001 020400 Census Tract 204, Autauga County, Alabama
5 01 001 020200 Census Tract 202, Autauga County, Alabama
6 01 001 020600 Census Tract 206, Autauga County, Alabama
P001001 H010001
1 1912 1912
2 10766 10585
3 3373 3373
4 4386 4386
5 2170 1989
6 3668 3668Notice that you had to specify the region, and since tracts are nested within states, you had to add regionin = stateget.
Two packages are available to access data gathered by the World Bank, {data360r} and {wbstats}. Install both packages, and then we can start with {data360r}, seeing what indicators, basic country-level information, and data-sets are available.
{data360r} seems moribund at the moment so the code shown below is not executed.
library(data360r)
#get all indicator metadata in Govdata360
get_metadata360(
site = "gov",
metadata_type = "indicators"
) -> df_indicators
#get all country metadata in TCdata360
get_metadata360(
metadata_type = 'countries'
) -> df_countries
#get all dataset metadata in TCdata360
get_metadata360(
metadata_type = 'datasets'
) -> df_datasets Once you have identified a particular table, note it’s ID, and then pull it. Say I want indicator 90 = Can a married woman register a business in the same way as a married man? I can do the same thing for more than one indicator by specifying the IDs.
get_data360(
indicator_id = c(90)
) -> df_ind90
get_data360(
indicator_id = c(28130, 28131)
) -> df_indtwo If I only want all data for a specific country or just the indicators we pull in df_ind1 for a specific country you could do:
get_data360(
country_iso3 = "IND"
) -> df_allone
search_360(
"woman business",
search_type = "indicator",
limit_results = 5
) -> df_ind1
get_data360(
indicator_id = df_ind1$id,
country_iso3 = "IND"
) -> df_allindtwo Now an example with two measures – legal age of marriage for boys and girls. Note that the package allows you to specify the long format (preferred) than the default wide format you see earlier results being returned in. Note also the use of timeframes = c() that allows you to specify the specific time period you want the indicators for.
search_360("marriage", search_type = "indicator")
get_data360(
indicator_id = c(204, 205),
timeframes = c(2016),
output_type = 'long'
) -> df_marriage
ggplot(
df_marriage,
aes(
x = Observation,
group = Indicator,
fill = Indicator
)
) +
geom_bar() +
theme(legend.position = "none") +
facet_wrap(~ Indicator, ncol = 1) +
labs(
x = "Legal Age of Marriage",
y = "Frequency",
title = "Legal Age of Marriage for Boys vs. Girls",
subtitle = "(2016)",
caption = "Source: World Bank Data"
)The {wbstats} package does pretty much the same thing. Let us see the core functionality by loading the library and then seeing what is available in terms of indicators, topics, and so on. We can then set the most current list of information in wb_cachelist to be used via new_cache. Doing so speeds up the operations and ensures that you are getting the latest data.
library(wbstats)
str(wb_cachelist, max.level = 1)List of 8
$ countries : tibble [297 × 18] (S3: tbl_df/tbl/data.frame)
$ indicators : tibble [16,643 × 8] (S3: tbl_df/tbl/data.frame)
$ sources : tibble [62 × 9] (S3: tbl_df/tbl/data.frame)
$ topics : tibble [21 × 3] (S3: tbl_df/tbl/data.frame)
$ regions : tibble [42 × 4] (S3: tbl_df/tbl/data.frame)
$ income_levels: tibble [7 × 3] (S3: tbl_df/tbl/data.frame)
$ lending_types: tibble [4 × 3] (S3: tbl_df/tbl/data.frame)
$ languages : tibble [23 × 3] (S3: tbl_df/tbl/data.frame)wb_cache() -> new_cache What indicators are available on the topic of corruption?
wb_search(pattern = "corruption") -> corruption_vars
head(corruption_vars)# A tibble: 6 × 3
indicator_id indicator indicator_desc
<chr> <chr> <chr>
1 CC.EST Control of Corruption: Estimate "Control of C…
2 CC.NO.SRC Control of Corruption: Number of So… "Control of C…
3 CC.PER.RNK Control of Corruption: Percentile R… "Control of C…
4 CC.PER.RNK.LOWER Control of Corruption: Percentile R… "Control of C…
5 CC.PER.RNK.UPPER Control of Corruption: Percentile R… "Control of C…
6 CC.STD.ERR Control of Corruption: Standard Err… "Control of C…If I want information from a particular source, say from Bloomberg,
wb_search(
pattern = "Bloomberg",
fields = "source_org"
) -> blmbrg_vars
head(blmbrg_vars)# A tibble: 2 × 3
indicator_id indicator indicator_desc
<chr> <chr> <chr>
1 GFDD.OM.02 Stock market return (%, year-on-year) Stock market ret…
2 GFDD.SM.01 Stock price volatility Stock price vola…Searching for indicators tied to multiple subjects is easy as well:
wb_search(
pattern = "poverty | unemployment | employment"
) -> povemply_vars
head(povemply_vars)# A tibble: 6 × 3
indicator_id indicator indicator_desc
<chr> <chr> <chr>
1 1.0.HCount.1.90usd Poverty Headcount ($1.90 a day) The poverty h…
2 1.0.HCount.2.5usd Poverty Headcount ($2.50 a day) The poverty h…
3 1.0.HCount.Mid10to50 Middle Class ($10-50 a day) Hea… The poverty h…
4 1.0.HCount.Ofcl Official Moderate Poverty Rate-… The poverty h…
5 1.0.HCount.Poor4uds Poverty Headcount ($4 a day) The poverty h…
6 1.0.HCount.Vul4to10 Vulnerable ($4-10 a day) Headco… The poverty h…Once we identify what we want, downloading the data is easy as well, needing us to specify just the indicator(s) and then the start and end dates, and then specific country codes if you want data for specific countries. Below I am pulling total population.
wb_data(
indicator = "SP.POP.TOTL",
start_date = 1960,
end_date = 2016
) -> pop_data1
wb_data(
country = c("ABW","AF", "SSF", "ECA", "IND", "CHN"),
indicator = "SP.POP.TOTL",
start_date = 1960,
end_date = 2016
) -> pop_data2
wb_data(
country = c("ABW","AF", "SSF", "ECA", "IND", "CHN"),
indicator = c("SP.POP.TOTL", "NY.GDP.MKTP.CD"),
start_date = 1960,
end_date = 2016
) -> pop_data3
head(pop_data3)# A tibble: 6 × 6
iso2c iso3c country date NY.GDP.MKTP.CD SP.POP.TOTL
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 AW ABW Aruba 1960 NA 54208
2 AW ABW Aruba 1961 NA 55434
3 AW ABW Aruba 1962 NA 56234
4 AW ABW Aruba 1963 NA 56699
5 AW ABW Aruba 1964 NA 57029
6 AW ABW Aruba 1965 NA 57357wb_data(
country = c("ABW","AF", "SSF", "ECA", "IND", "CHN"),
indicator = c("SP.POP.TOTL", "NY.GDP.MKTP.CD"),
start_date = 1960,
end_date = 2016,
return_wide = TRUE
) -> pop_data4
head(pop_data4)# A tibble: 6 × 6
iso2c iso3c country date NY.GDP.MKTP.CD SP.POP.TOTL
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 AW ABW Aruba 1960 NA 54208
2 AW ABW Aruba 1961 NA 55434
3 AW ABW Aruba 1962 NA 56234
4 AW ABW Aruba 1963 NA 56699
5 AW ABW Aruba 1964 NA 57029
6 AW ABW Aruba 1965 NA 57357By default wb_data() will return the data in long format but not necessarily in a tidy format. If you want the data returned on call in a wide format, specify return_wide = TRUE and you will have tidy data.
If you will be working with dates, whether for plotting purposes or otherwise, then activate the date_as_class_date = TRUE switch. Otherwise you will have to manually format the date variables.
wb_data(
country = c("IND", "CHN"),
indicator = c("SP.POP.TOTL"),
start_date = 1960,
end_date = 2016,
return_wide = TRUE,
date_as_class_date = TRUE
) -> pop_data5
head(pop_data5)# A tibble: 6 × 10
iso2c iso3c country date SP.POP.TOTL unit obs_status footnote
<chr> <chr> <chr> <date> <dbl> <chr> <chr> <chr>
1 CN CHN China 1960-01-01 667070000 <NA> <NA> <NA>
2 CN CHN China 1961-01-01 660330000 <NA> <NA> <NA>
3 CN CHN China 1962-01-01 665770000 <NA> <NA> <NA>
4 CN CHN China 1963-01-01 682335000 <NA> <NA> <NA>
5 CN CHN China 1964-01-01 698355000 <NA> <NA> <NA>
6 CN CHN China 1965-01-01 715185000 <NA> <NA> <NA>
# … with 2 more variables: last_updated <date>, obs_resolution <chr>library(scales)
ggplot(
pop_data5,
aes(
x = date,
y = SP.POP.TOTL)
) +
geom_line(
aes(
color = country
)
) +
scale_y_continuous(labels = comma) +
scale_x_date(date_breaks = "10 years") +
theme(legend.position = "bottom") +
labs(
x = "Date",
y = "Total Population"
) +
theme_minimal()In the case that you want the most recent data, use mrv switch, and note that the number indicates how many of the most recent values you would like to access. This will also return rows of data where there are no recent values. If you wish to avoid this latter result, use mrnev instead.
wb_data(
country = "all",
indicator = c("SP.POP.TOTL"),
mrv = 1,
return_wide = TRUE,
date_as_class_date = TRUE
) -> pop_data6
head(pop_data6)# A tibble: 6 × 10
iso2c iso3c country date SP.POP.TOTL unit obs_status footnote
<chr> <chr> <chr> <date> <dbl> <chr> <chr> <chr>
1 AW ABW Aruba 2020-01-01 106766 <NA> <NA> <NA>
2 <NA> AFE Africa… 2020-01-01 677243299 <NA> <NA> <NA>
3 AF AFG Afghan… 2020-01-01 38928341 <NA> <NA> <NA>
4 <NA> AFW Africa… 2020-01-01 458803476 <NA> <NA> <NA>
5 AO AGO Angola 2020-01-01 32866268 <NA> <NA> <NA>
6 AL ALB Albania 2020-01-01 2837743 <NA> <NA> Extrapo…
# … with 2 more variables: last_updated <date>, obs_resolution <chr>Eritrea has missing data. To avoid this we could have done the following (code not executed below) …
wb_data(
country = "all",
indicator = c("SP.POP.TOTL"),
mrnev = 1,
freq = "Y",
return_wide = TRUE,
date_as_class_date = TRUE
) -> pop_data7
head(pop_data7)We could have also asked for data by countries only, administrative regions, and so on.
wb_data(
country = "regions_only",
indicator = c("SP.POP.TOTL"),
mrv = 1
) -> regions
wb_data(
country = "admin_regions_only",
indicator = c("SP.POP.TOTL"),
mrv = 1
) -> admin_regions
wb_data(
country = "income_levels_only",
indicator = c("SP.POP.TOTL"),
mrv = 1
) -> income_levels
wb_data(
country = "lending_types_only",
indicator = c("SP.POP.TOTL"),
mrv = 1
) -> lending_typesWhat if I want the data in Chinese?
wb_data(
country = "lending_types_only",
indicator = c("SP.POP.TOTL"),
mrv = 1,
lang = "zh"
) -> lending_types_chinese
head(lending_types_chinese)# A tibble: 3 × 9
iso2c iso3c country date SP.POP.TOTL unit obs_status footnote
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
1 XF IBD 只有IBRD 2020 4862388283 <NA> <NA> <NA>
2 XH IDB IDA混合 2020 574159138 <NA> <NA> <NA>
3 XI IDX 只有IDA 2020 1134444535 <NA> <NA> <NA>
# … with 1 more variable: last_updated <date>The {dataRetrieval} package gives you easy access to water data gathered and warehoused by the USGS, USDA, EPA, and other entities. The package has an excellent tutorial available here so I will not go into too many details and nuances here. Start by installing the package and then loading it.
library(dataRetrieval)You will need to know the site(s) you are interested in, as well as the parameters and statistics of interest. The package comes with a built-in data-set that shows you the parameters available, and the complete list of statistics is available here. Sites can be located from this inventory.
parameterCdFile -> parameterCdFile
names(parameterCdFile)[1] "parameter_cd" "parameter_group_nm" "parameter_nm"
[4] "casrn" "srsname" "parameter_units" If you are curious about a specific parameter, you can see what all is available for it. I’ll look for anything related to the keyword storm, and also what parameter units are available.
parameterCdFile[
grep("storm",
parameterCdFile$parameter_nm,
ignore.case = TRUE),
] -> stormq
unique(stormq$parameter_units)[1] "nu" "hours" "minutes" "Mgal" "ft3/s" "mgd"
[7] "in" Let us do a quick grab of some data for the Hocking River at Athens, Ohio.
siteNo <- "03159500"
pCode <- "00065"
start.date <- "2014-10-01"
end.date <- "2021-03-20"
readNWISuv(
siteNumbers = siteNo,
parameterCd = pCode,
startDate = start.date,
endDate = end.date
) -> hocking Since the column names are based on parameter codes and hence cryptic, you can clean them up, and also see other attributes embedded in the data-set.
names(hocking)[1] "agency_cd" "site_no" "dateTime"
[4] "X_00065_00000" "X_00065_00000_cd" "tz_cd" renameNWISColumns(hocking) -> hocking
names(hocking)[1] "agency_cd" "site_no" "dateTime" "GH_Inst" "GH_Inst_cd"
[6] "tz_cd" names(attributes(hocking))[1] "names" "row.names" "class" "url"
[5] "siteInfo" "variableInfo" "disclaimer" "statisticInfo"
[9] "queryTime" If I wanted data for multiple sites, I could find the site numbers and then grab the data.
sites <- c("03158200", "03159246")
pcode <- "00065"
start.date <- "2014-10-01"
end.date <- "2025-03-21"
readNWISuv(
siteNumbers = sites,
parameterCd = pcode,
startDate = start.date,
endDate = end.date
) -> hocking2
renameNWISColumns(hocking2) -> hocking2 Now a simple time-series plot of gage height for both sites. Note that although I asked for data going far back, not all sites have all data for all time periods so it helps to check the site inventory first.
attr(hocking2, "variableInfo") -> parameterInfo
ifelse(
hocking2$site_no == "03158200",
"Monday Creek at Doanville",
"Sunday Creek Below Millfield"
) -> hocking2$station
as.Date(
as.character(hocking2$dateTime),
format = "%Y-%m-%d"
) -> hocking2$mydates
ggplot(
data = hocking2,
aes(
x = mydates,
y = GH_Inst,
color = station
)
) +
geom_line() +
labs(
x = "",
y = parameterInfo$variableDescription
) +
scale_x_date(date_breaks = "16 weeks") +
theme_minimal() +
theme(legend.position = "bottom") For attribution, please cite this work as
Ruhil (2025, March 23). Using APIs from the Census Bureau & other sources. Retrieved from https://aniruhil.org/courses/mpa6020/handouts/module06.html
BibTeX citation
@misc{ruhil2025using,
author = {Ruhil, Ani},
title = {Using APIs from the Census Bureau & other sources},
url = {https://aniruhil.org/courses/mpa6020/handouts/module06.html},
year = {2025}
}