2.9 Buildings

2018 data from the local County Assessor’s Office includes property use type and total building square footage for every parcel in Stockton. Because the assessor’s record of total building square footage is of poor quality, we chose to instead use Microsoft’s U.S. Building Footprint Data, joined each building footprint to a parcel (thereby also assigning a residential or commercial type based on the zoning of the parcel), and multiplied this footprint area by number of stories, which came either from the assessor’s record or from the Microsoft data. This dataset may still contain significant errors (e.g. mixed use parcels are mis-identified, the stories record may be incorrect, Microsoft’s building footprints may be incorrect, etc.), but we decided it was the best representation of total building square footage in Stockton.

The following maps show the density of residents per acre across Stockton, followed by the density of jobs per acre of commercials buildings, at the block group level. Residential population was obtained from ACS 2017 5-yr data, and job count was obtained from LODES 2017.

pge_stockton_filtered <- 
  pge_data %>% 
  filter(!CUSTOMERCLASS %in% c("EI","EA") & (YEAR == 2018)) %>%
  dplyr::select(-YEAR,-TOTALTCO2E,-TCO2EPERCUST)

zips_spread <- 
  pge_stockton_filtered %>% 
  group_by(ZIPCODE, CUSTOMERCLASS) %>% 
  summarise(
    TOTALKBTU = sum(TOTALKBTU, na.rm=T),
    TOTALCUSTOMERS = mean(TOTALCUSTOMERS, na.rm=T)
  ) %>% 
  mutate(KBTUPERCUST = TOTALKBTU/TOTALCUSTOMERS) %>% 
  gather(key = "type", value, TOTALKBTU:KBTUPERCUST) %>% 
  unite(temp,CUSTOMERCLASS,type) %>% 
  spread(temp,value) 

#next, since the previous line fully disaggregates by type AND class, but i also consider it valuable to disaggregate by type OR class individually on their own, i do a second and third spread. all of these "spreads" will be joined back to an original in the final step.
zips_spread_2 <- 
  pge_stockton_filtered %>% 
  group_by(ZIPCODE, CUSTOMERCLASS) %>% 
  summarise(
    TOTALKBTU = sum(TOTALKBTU, na.rm=T),
    TOTALCUSTOMERS = mean(TOTALCUSTOMERS, na.rm=T)
  ) %>% 
  ungroup() %>% 
  mutate(ENERGYTYPE = substr(CUSTOMERCLASS,1,1)) %>% 
  dplyr::select(-CUSTOMERCLASS) %>% 
  group_by(ZIPCODE, ENERGYTYPE) %>% 
  summarise(
    TOTALKBTU = sum(TOTALKBTU, na.rm=T),
    TOTALCUSTOMERS = sum(TOTALCUSTOMERS, na.rm=T)
  ) %>% 
  mutate(KBTUPERCUST = TOTALKBTU/TOTALCUSTOMERS) %>% 
  gather(key = "type", value, TOTALKBTU:KBTUPERCUST) %>% 
  unite(temp,ENERGYTYPE,type) %>% 
  spread(temp,value)

#note that when combining into SECTOR, say electricity and gas for commercial, the summary function for TOTALCUSTOMERS switches to max(), with the reasoning that many of these commercial customers have both electricity and gas, so we would assume that the correct total # of customers is closer to the max of either, than to add them together (which presumes that there are no overlapping electricity and gas customers).
zips_spread_3 <- 
  pge_stockton_filtered %>%
  group_by(ZIPCODE, CUSTOMERCLASS) %>% 
  summarise(
    TOTALKBTU = sum(TOTALKBTU, na.rm=T),
    TOTALCUSTOMERS = mean(TOTALCUSTOMERS, na.rm=T)
  ) %>% 
  ungroup() %>% 
  mutate(SECTOR = substr(CUSTOMERCLASS,2,2)) %>% 
  dplyr::select(-CUSTOMERCLASS) %>% 
  group_by(ZIPCODE, SECTOR) %>% 
  summarise(TOTALKBTU = sum(TOTALKBTU, na.rm=T),
            TOTALCUSTOMERS = max(TOTALCUSTOMERS, na.rm=T)) %>% 
  mutate(KBTUPERCUST = TOTALKBTU/TOTALCUSTOMERS) %>% 
  gather(key = "type", value, TOTALKBTU:KBTUPERCUST) %>% 
  unite(temp,SECTOR,type) %>% 
  spread(temp,value)

#the final step here is to get the actual total totals for kbtu and TCO2E indicators, and then join all the "subtotals" from the previous 3 spreads.
zips_total <- 
  pge_stockton_filtered %>% 
  group_by(ZIPCODE) %>% 
  summarize(TOTALKBTU = sum(TOTALKBTU, na.rm=T)) %>% 
  left_join(zips_spread_3, by = "ZIPCODE") %>% 
  mutate(TOTALCUSTOMERS = R_TOTALCUSTOMERS + C_TOTALCUSTOMERS) %>% 
  left_join(zips_spread_2, by = "ZIPCODE") %>% 
  left_join(zips_spread, by = "ZIPCODE")

#join this massive summary back to the zcta geometries
zcta_stockton_joined <- 
  zcta_stockton %>% 
  left_join(zips_total, by="ZIPCODE")

#below I've brought in all the code from stockton_bldg.R, but it can all be skipped by going to the load() of stockton_bldg.Rdata at the bottom. Kevin, whatever refinements you've made, go ahead and make them here.

# sjc_bldg <- 
#   read_csv("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/USBuildingFootprints/ca_06077_footprints.csv") %>% 
#   st_as_sf(wkt = "WKT") %>% 
#   st_set_crs(4326) %>% 
#   st_transform(st_crs(zcta_stockton)) %>% 
#   mutate(id = row_number())
# 
# stockton_bldg <- 
#   sjc_bldg[which(sjc_bldg$id %in% st_centroid(sjc_bldg)[stockton_boundary_influence,]$id),]
# 
# sjc_parcels <- 
#   st_read("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/Parcels/Parcels.shp") %>% 
#   st_transform(st_crs(zcta_stockton)) %>% 
#   st_make_valid()
# 
# save(sjc_parcels, file = "C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/Parcels/sjc_parcels.Rdata")
# 
# stockton_parcels <- 
#   sjc_parcels[stockton_boundary_influence,]
# 
# bldg_parcel_join <- 
#   st_join(st_centroid(stockton_bldg), stockton_parcels) %>%
#   dplyr::select(APN, id, STAREA__) %>% 
#   rename(area = STAREA__) %>% 
#   st_set_geometry(NULL)
# 
# sjc_zoning <- 
#   st_read("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/Zoning/Zoning.shp") %>% 
#   st_transform(st_crs(zcta_stockton)) %>% 
#   filter(ZNLABEL != "STOCKTON")
# 
# stockton_zoning <- 
#   st_read("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/Stockton_Zoning/Zoning.shp") %>% 
#   st_transform(st_crs(zcta_stockton))
# 
# bldg_zoning_join <- 
#   st_join(st_centroid(stockton_bldg), stockton_zoning) %>% 
#   dplyr::select(ZONE, id) %>% 
#   st_set_geometry(NULL)
# 
# bldg_zoning_join_uninc <- 
#   st_join(st_centroid(stockton_bldg), sjc_zoning) %>% 
#   dplyr::select(ZNCODE, id) %>% 
#   st_set_geometry(NULL)
# 
# bldg_zoning_join %<>% 
#   merge(bldg_zoning_join_uninc) %>% 
#   mutate(ZONE = ifelse(is.na(ZONE),as.character(ZNCODE),as.character(ZONE))) %>% 
#   dplyr::select(ZONE, id)
# 
# sjc_bgs <- 
#   block_groups("California", "San Joaquin County")
# 
# stockton_bgs <- 
#   sjc_bgs[stockton_boundary_influence,]
# 
# bldg_bg_join <- 
#   st_join(st_centroid(stockton_bldg), stockton_bgs) %>% 
#   dplyr::select(GEOID, id) %>% 
#   st_set_geometry(NULL)
# 
# bldg_zcta_join <- 
#   st_join(st_centroid(stockton_bldg), zcta_stockton) %>% 
#   dplyr::select(ZIPCODE, id) %>% 
#   st_set_geometry(NULL)
# 
# load("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/tax_sjc.Rdata")
# 
# tax_sjc <- 
#   tax_sjc %>% 
#   mutate(APN = as.numeric(apn)) %>% 
#   st_set_geometry(NULL)
# 
# stockton_bldg_final <-
#   stockton_bldg %>%
#   left_join(bldg_parcel_join, by="id") %>%
#   left_join(bldg_zoning_join, by="id") %>%
#   left_join(bldg_bg_join, by="id") %>%
#   left_join(bldg_zcta_join, by="id") %>%
#   left_join(tax_sjc, by="APN") %>%
#   mutate( # the following takes the bldg data, computes bldgground (ground footprint, where the factor conversion is sqm to sqft) and bldgtot (total sqft).
#     bldgground = as.numeric(st_area(WKT))*10.7639,
#     bldgtot =
#       ifelse(
#         is.na(stories),
#         bldgground,
#         bldgground*stories
#       ),
#     ZIPCODE = as.numeric(ZIPCODE),
#     ZONING = 
#       case_when(
#         substr(ZONE,1,1) == "R" ~ "R",
#         substr(ZONE,1,2) == "(R" ~ "R",
#         substr(ZONE,1,1) %in% c("C","P","M","U") ~ "C", 
#         substr(ZONE,1,2) %in% c("(C","(P") ~ "C",
#         TRUE ~ "O"
#       )
#   )
# 
# save(stockton_bldg_final, file = "C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/stockton_bldg.Rdata")
load("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/stockton_bldg.Rdata")
# 
# stockton_bldg_17 <- 
#   readOGR(dsn = "C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/USBuildingFootprints/Stockton2017", layer = "Stockton") %>% 
#   st_as_sf() %>% 
#   st_transform(st_crs(stockton_bldg_final))
# 
# save(stockton_bldg_17, file = "C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/USBuildingFootprints/stockton_bldg_17.Rdata")
# load("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/USBuildingFootprints/stockton_bldg_17.Rdata")

# # match height from 2017 data
# bldg_join <-
#   stockton_bldg_final %>%
#   st_join(
#     stockton_bldg_17 %>%
#       st_transform(26910) %>%
#       st_buffer(-1) %>%
#       st_transform(st_crs(stockton_bldg_final))
#   ) %>%
#   group_by(id) %>%
#   summarize(
#     Height = mean(Height,na.rm=T)
#   ) %>%
#   left_join(stockton_bldg_final %>% st_set_geometry(NULL),by="id")
# 
# bldg_join <-
#   bldg_join %>% 
#   mutate(
#     stories_by_height =
#       case_when(
#         GEOID == "060770039001" ~ 1, #Naval base
#         is.na(Height) & is.na(stories) ~ 1,
#         is.na(Height) ~ round(stories),
#         !is.na(stories) ~ round(stories),
#         Height < mean(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==1),]$Height) + sd(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==1),]$Height) ~ 1,
#         Height < mean(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==2),]$Height) + sd(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==2),]$Height) ~ 2,
#         Height < mean(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==3),]$Height) + sd(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==3),]$Height) ~ 3,
#         Height < mean(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==4),]$Height) + sd(bldg_join[which(!is.na(bldg_join$Height)&round(bldg_join$stories)==4),]$Height) ~ 4,
#         TRUE ~ round(Height/5)
#       ),
#     bldg_tot =
#       bldgground*stories_by_height
#   )
# # 6033 unmatched when using -1m buffer. #2637 of those did not have assessor stories data, set to 1F. Of other 3396, 2090 1F, 1248 2F, 58 3F.
# # 70757 had both 2017 heights and assessor stories data, go with assessor stories data. 51676 had assessor stories = 1, kept at 1F. 18974 had assessor stories = 2, kept at 2F. 107 had assessor stories > 2, set to 2.5F.
# # mean height for assessor 1F: 11.30, sd 1.63. mean height for assessor 2F: 13.98, sd 2.25. mean height for assessor 3F: 15.39, sd 2.99. mean height for assessor 4F: 14.42, sd 5.35.
# # 24660 had no assessor stories data but had 2017 heights, use mean + 1 sd from above as thresholds. 1F: 12.92. 2F: 16.23. 3F: 18.38. 4F: 19.77
# # 15991 1F, 5642 2F, 908 3F, 1512 4F, 541 5F, 29 6F, 12 7F, 15 8F, 4 9F, 5 10F, 1 12F.
# 
# save(bldg_join, file = "C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/USBuildingFootprints/bldg_join.Rdata")
load("C:/Users/derek/Google Drive/City Systems/Stockton Green Economy/USBuildingFootprints/bldg_join.Rdata")

bg_bldg <-
  bldg_join %>% 
  st_set_geometry(NULL) %>% 
  group_by(GEOID, ZONING) %>% 
  summarize(
    numbldg = n(),
    bldgtot = sum(bldgtot, na.rm=T)
  ) %>% 
  ungroup() %>% 
  left_join(stockton_bgs_full, by = "GEOID") %>% 
  st_as_sf() %>% 
  filter(!is.na(st_dimension(.)))

#res/sqft per block group
bg_pop <-
  getCensus(
    name = "acs/acs5",
    vintage = 2017, # satellite imagery is from around 2014/2015
    vars = c("B01003_001E"),
    region = "block group:*",
    regionin = "state:06+county:077"
  ) %>%
  transmute(
    GEOID = paste0(state,county,tract,block_group),
    pop = B01003_001E
  )

bg_bldg_norm_r <-
  bg_bldg %>% 
  filter(ZONING == "R") %>% 
  left_join(bg_pop, by = "GEOID") %>% 
  mutate(
    `Residents per 1,000 sqft` = pop/bldgtot*1000,
    residential_density = bldgtot/(as.numeric(st_area(.))/4046.86)
  )

# map = mapview(bg_bldg_norm_r[,c("residential_density")],zcol=c("residential_density"), legend = TRUE, layer.name = "Residential</br>density</br>per acre")
# 
# mapshot(map, url = "map-residential-density.html")

knitr::include_url("https://city.systems/stockton-greeneconomy/map-residential-density.html")

Figure 2.38: Stockton residents per 1,000 sqft of residential building area, per block group.

\(~\)

# bldg_norm_r <-
#   bg_bldg %>% 
#   filter(ZONING == "R") %>% 
#   left_join(bg_pop, by = "GEOID") %>% 
#   summarize(
#     pop = sum(pop, na.rm=T),
#     bldgtot = sum(bldgtot, na.rm=T)
#   ) %>% 
#   mutate(`Residents/1000sqft` = pop/bldgtot*1000)
# 
# bldg_norm_c <-
#   bg_bldg %>% 
#   filter(ZONING == "C") %>% 
#   left_join(
#     stockton_wac %>% 
#       filter(year == 2017) %>% 
#       dplyr::select(jobs = C000, GEOID), 
#     by = "GEOID"
#   ) %>% 
#   summarize(
#     jobs = sum(jobs, na.rm=T),
#     bldgtot = sum(bldgtot, na.rm=T)
#   ) %>% 
#   mutate(`Jobs/1000sqft` = jobs/bldgtot*1000)

bg_bldg_norm_c <-
  bg_bldg %>% 
  filter(ZONING == "C") %>% 
  left_join(
    stockton_wac %>% 
      as.data.frame() %>% 
      filter(year == 2015) %>% 
      dplyr::select(jobs=C000,GEOID), 
    by = "GEOID"
  ) %>% 
  mutate(
    `Jobs per 1,000 sqft` = jobs/bldgtot*1000,
    commercial_density = bldgtot/(as.numeric(st_area(.))/4046.86)
  ) %>% 
  filter(!GEOID %in% c("060770039001","060770039002"))

# map = mapview(bg_bldg_norm_c[,c("commercial_density")],zcol=c("commercial_density"), legend = TRUE, layer.name = "Commercial</br>density</br>per acre")
# 
# mapshot(map, url = "map-commercial-density.html")

knitr::include_url("https://city.systems/stockton-greeneconomy/map-commercial-density.html")

Figure 2.39: Stockton jobs per 1,000 sqft of commercial building area, per block group.

\(~\)

Because we do not have robust historical building datasets, we cannot analyze the change in building square footages over time. For forecasting, we assumed an overall fixed building utilization estimate of 1.7 residents per 1000 sqft of residential buildings, and 1.3 jobs per 1000 sqft of commercial buildings. In the intervention stage, if we expect a strategy would change building utilization, we would make a reduction in energy consumption relative to these baseline utilization rates.

As one last couple of building-specific analyses, we aggregated building square footage information up to the zipcode level and compared with energy usage data to see residential and commercial building energy efficiency at a recent snapshot of time, spatially distributed across Stockton.

#disaggregating by zipcode and zoning type of the buildings to get summaries of # of buildings, total ground footprint, and total bldg sqft. the filter removes NAs in the data, which could be missing zones in specific zipcodes.
zcta_bldg_spread <- 
  stockton_bldg_final %>% 
  st_set_geometry(NULL) %>% 
  group_by(ZIPCODE, ZONING) %>% 
  summarise(
    NUMBLDG = n(),
    TOTSQFTGROUND = sum(bldgground, na.rm=T),
    TOTSQFT = sum(bldgtot, na.rm=T)
  ) %>% 
  filter(!is.na(ZIPCODE) & !is.na(ZONING)) %>% 
  gather(key, value, NUMBLDG:TOTSQFT) %>% 
  unite(temp,ZONING,key) %>% 
  spread(temp,value)

#summarize totals by zipcode, then join the zoning-specific subtotals to it.
zcta_bldg_summary <-
  stockton_bldg_final %>% 
  st_set_geometry(NULL) %>% 
  group_by(ZIPCODE) %>% 
  summarise(
    NUMBLDG = n(),
    TOTSQFTGROUND = sum(bldgground, na.rm=T),
    TOTSQFT = sum(bldgtot, na.rm=T)
  ) %>% 
  left_join(zcta_bldg_spread, by = "ZIPCODE")

zcta_bldg_stockton_joined <- 
  zcta_stockton_joined %>% 
  filter(ZIPCODE != 95211) %>% 
  left_join(zcta_bldg_summary, by = "ZIPCODE") %>% 
  mutate(
    ER_KBTUperSQFT = ER_TOTALKBTU/R_TOTSQFT, 
    GR_KBTUperSQFT = GR_TOTALKBTU/R_TOTSQFT, 
    EC_KBTUperSQFT = EC_TOTALKBTU/C_TOTSQFT, 
    GC_KBTUperSQFT = GC_TOTALKBTU/C_TOTSQFT, 
    R_KBTUperSQFT = R_TOTALKBTU/R_TOTSQFT, 
    C_KBTUperSQFT = C_TOTALKBTU/C_TOTSQFT, 
    E_KBTUperSQFT = E_TOTALKBTU/TOTSQFT, 
    G_KBTUperSQFT = G_TOTALKBTU/TOTSQFT, 
    KBTUperSQFT = TOTALKBTU/TOTSQFT
  )

# map = mapview(zcta_bldg_stockton_joined[,c("R_KBTUPERCUST")], zcol=c("R_KBTUPERCUST"), legend = TRUE, layer.name = "kBTU per customer,</br>Residential")
# 
# mapshot(map, url = "map-res-customer.html")

knitr::include_url("https://city.systems/stockton-greeneconomy/map-res-customer.html")

Figure 2.40: Stockton residential kBTU per sqft of residential building area, per zip code.

\(~\)

# map = mapview(zcta_bldg_stockton_joined[,c("C_KBTUPERCUST")], zcol=c("C_KBTUPERCUST"), legend = TRUE, layer.name = "kBTU per customer,</br>Commercial")
# 
# mapshot(map, url = "map-comm-customer.html")

knitr::include_url("https://city.systems/stockton-greeneconomy/map-comm-customer.html")

Figure 2.41: Stockton commercial kBTU per sqft of commercial building area, per zip code.

\(~\)

zcta_bldg_stockton_summary <- 
  zcta_bldg_stockton_joined %>% 
  st_set_geometry(NULL) %>% 
  summarise_at(
    c(
      "NUMBLDG", 
      "R_NUMBLDG", 
      "C_NUMBLDG",
      "TOTSQFTGROUND",
      "TOTSQFT",
      "R_TOTSQFT",
      "C_TOTSQFT",
      "R_TOTSQFTGROUND",
      "C_TOTSQFTGROUND",
      "TOTALKBTU",
      "E_TOTALKBTU",
      "G_TOTALKBTU",
      "R_TOTALKBTU",
      "C_TOTALKBTU", 
      "ER_TOTALKBTU", 
      "GR_TOTALKBTU", 
      "EC_TOTALKBTU", 
      "GC_TOTALKBTU"
    ), sum, na.rm=T
  ) %>% 
  mutate(
    ER_KBTUperSQFT = ER_TOTALKBTU/R_TOTSQFT, 
    GR_KBTUperSQFT = GR_TOTALKBTU/R_TOTSQFT, 
    EC_KBTUperSQFT = EC_TOTALKBTU/C_TOTSQFT, 
    GC_KBTUperSQFT = GC_TOTALKBTU/C_TOTSQFT, 
    R_KBTUperSQFT = R_TOTALKBTU/R_TOTSQFT, 
    C_KBTUperSQFT = C_TOTALKBTU/C_TOTSQFT, 
    E_KBTUperSQFT = E_TOTALKBTU/TOTSQFT, 
    G_KBTUperSQFT = G_TOTALKBTU/TOTSQFT, 
    KBTUperSQFT = TOTALKBTU/TOTSQFT,
  )

zcta_bldg_stockton_table <-
  data.frame(
    Type = c(
      "Residential Electricity",
      "Residential Gas",
      "Commercial Electricity",
      "Commercial Gas",
      "Residential Total",
      "Commercial Total",
      "Electricity Total",
      "Gas Total",
      "Total"
    ),
    Bldg = c(
      NA,
      NA,
      NA,
      NA,
      zcta_bldg_stockton_summary$R_NUMBLDG,
      zcta_bldg_stockton_summary$C_NUMBLDG,
      NA,
      NA,
      zcta_bldg_stockton_summary$R_NUMBLDG+zcta_bldg_stockton_summary$C_NUMBLDG
    ),
    Footprint = c(
      NA,
      NA,
      NA,
      NA,
      zcta_bldg_stockton_summary$R_TOTSQFTGROUND,
      zcta_bldg_stockton_summary$C_TOTSQFTGROUND,
      NA,
      NA,
      zcta_bldg_stockton_summary$R_TOTSQFTGROUND+zcta_bldg_stockton_summary$C_TOTSQFTGROUND
    ),
    sqft = c(
      NA,
      NA,
      NA,
      NA,
      zcta_bldg_stockton_summary$R_TOTSQFT,
      zcta_bldg_stockton_summary$C_TOTSQFT,
      NA,
      NA,
      zcta_bldg_stockton_summary$R_TOTSQFT+zcta_bldg_stockton_summary$C_TOTSQFT
    ),
    kbtu = c(
      zcta_bldg_stockton_summary$ER_TOTALKBTU,
      zcta_bldg_stockton_summary$GR_TOTALKBTU,
      zcta_bldg_stockton_summary$EC_TOTALKBTU,
      zcta_bldg_stockton_summary$GC_TOTALKBTU,
      zcta_bldg_stockton_summary$R_TOTALKBTU,
      zcta_bldg_stockton_summary$C_TOTALKBTU,
      zcta_bldg_stockton_summary$E_TOTALKBTU,
      zcta_bldg_stockton_summary$G_TOTALKBTU,
      zcta_bldg_stockton_summary$TOTALKBTU
    ),
    kbtusqft = c(
      zcta_bldg_stockton_summary$ER_KBTUperSQFT,
      zcta_bldg_stockton_summary$GR_KBTUperSQFT,
      zcta_bldg_stockton_summary$EC_KBTUperSQFT,
      zcta_bldg_stockton_summary$GC_KBTUperSQFT,
      zcta_bldg_stockton_summary$R_KBTUperSQFT,
      zcta_bldg_stockton_summary$C_KBTUperSQFT,
      zcta_bldg_stockton_summary$E_KBTUperSQFT,
      zcta_bldg_stockton_summary$G_KBTUperSQFT,
      zcta_bldg_stockton_summary$TOTALKBTU/(zcta_bldg_stockton_summary$R_TOTSQFT+zcta_bldg_stockton_summary$C_TOTSQFT)
    )
  ) %>% 
  transmute(
    Type = Type,
    `Number of Buildings` = 
      ifelse(
        is.na(Bldg),
        "",
        prettyNum(round(Bldg,-3),big.mark = ",")
      ),
    `Total Footprint` =
      ifelse(
        is.na(Footprint),
        "",
        prettyNum(round(Footprint,-6),big.mark = ",")
      ),
    `Total Sqft` =
      ifelse(
        is.na(sqft),
        "",
        prettyNum(round(sqft,-6),big.mark = ",")
      ),
    `Total kBTU` = prettyNum(round(kbtu,-7),big.mark = ","),
    `kBTU/sqft` = prettyNum(round(kbtusqft),big.mark = ",")
  )

kable(
  zcta_bldg_stockton_table, 
  booktabs = TRUE,
  caption = 'Summary of Building count, size, and energy usage for residential/commercial and electricity/gas. Data from PG&E, 2018.'
) %>% 
  kable_styling() %>% 
  scroll_box(width = "100%")
Table 2.25: Summary of Building count, size, and energy usage for residential/commercial and electricity/gas. Data from PG&E, 2018.
Type Number of Buildings Total Footprint Total Sqft Total kBTU kBTU/sqft
Residential Electricity 2,700,000,000 9
Residential Gas 4,370,000,000 15
Commercial Electricity 2,620,000,000 45
Commercial Gas 1,210,000,000 21
Residential Total 91,000 237,000,000 288,000,000 7,080,000,000 25
Commercial Total 6,000 55,000,000 58,000,000 3,840,000,000 66
Electricity Total 5,330,000,000 13
Gas Total 5,580,000,000 14
Total 97,000 291,000,000 346,000,000 10,910,000,000 32

\(~\)