-
Notifications
You must be signed in to change notification settings - Fork 5
/
transform_frs.R
64 lines (51 loc) · 2.09 KB
/
transform_frs.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
# Transform EPA facility registry service data ---------------------------------
library(fs)
library(sf)
library(tidyverse)
# source functions
dir_ls(here::here("src/functions")) %>% walk(~source(.x))
# path to save raw data and standard projection
data_path <- Sys.getenv("WSB_DATA_PATH")
staging_path <- Sys.getenv("WSB_STAGING_PATH")
epsg <- as.numeric(Sys.getenv("WSB_EPSG"))
# Read un-zipped geodatabase (~6GB so querying on water to reduce file size)
# First look at available layers
frs_layers <- st_layers(dsn = path(data_path, "frs/FRS_INTERESTS.gdb"))
# SQL query to target facilities with water focus
get_water_frs <- "
SELECT *
FROM FACILITY_INTERESTS
WHERE INTEREST_TYPE IN (
'COMMUNITY WATER SYSTEM',
'NON-TRANSIENT NON-COMMUNITY WATER SYSTEM',
'TRANSIENT NON-COMMUNITY WATER SYSTEM',
'WATER TREATMENT PLANT',
'DRINKING WATER PROGRAM',
'DRINKING WATER SYSTEM'
)"
# Read layer for FRS_INTERESTS with conditional query on `INTEREST_TYPE`.
# Then, transform to standard epsg.
frs_water <- path(data_path, "frs/FRS_INTERESTS.gdb") %>%
st_read(query = get_water_frs,
layer = "FACILITY_INTERESTS",
stringsAsFactors = FALSE) %>%
st_transform(epsg)
cat("Read labeled FRS layer and transformed to CRS:", epsg, "\n ")
# Visualize points
#plot(st_geometry(frs_water), pch = 1, col = 'blue')
# General cleaning --------------------------------------------------------
# Set column names to lower case, clean names, clean whitespace,
# split PWSID and Facility ID from pgm_sys_id, add reported state name
frs_water <- frs_water %>%
rename(geometry = Shape) %>%
janitor::clean_names() %>%
f_clean_whitespace_nas() %>%
mutate(pwsid = word(pgm_sys_id, 1),
state = substr(pwsid, 1, 2),
facility_id = word(pgm_sys_id, 2),
facility_id = ifelse(pwsid == facility_id, NA, facility_id))
# Write to geopackage ----------------------------------------------------
path_out <- path(staging_path, "frs.gpkg")
if(file_exists(path_out)) file_delete(path_out)
st_write(frs_water, path_out)
cat("Wrote FRS data to geopackage. \n")