Sfs domain clipped rural boundary

sfs_domain_geojsons/define_SFS_domain_shapes.R

@@ -19,12 +19,23 @@ citytracts <- masterSpatialDB(shape_level = 'census_tract', source='seattle_geoj
 wa <- masterSpatialDB(shape_level = 'puma', source='wa_geojson')
 
 # surround pumas
-surround <- wa %>% filter(!(PUMACE10 %in% city$residence_puma))
+surround <- wa %>% filter(!(PUMACE10 %in% city$PUMA5CE))
+
+## merge neighborhoods in city with pumas in surround
+# keep state, name, lowest geoid, puma, tract, domain
+domain <- st_sf( regional_name = c(paste('Seattle--',as.character(city$NEIGHBO),sep=''), as.character(surround$NAME10)),
+                 STATE = 53,
+                 GEOID = c(rep('NA',nrow(city)), paste('53',surround$residence_puma,sep='')),
+                 domain = 'SFS_year2',
+                 geometry = c(city$geometry,surround$geometry)
+               )
 
 # zip maps
+zipsAll <- zctas(cb=TRUE, starts_with = '98')
+zipsAll <- st_transform(st_as_sf(zipsAll),4326)
+
 dat<- read_excel('Zip_Codes_SFS-Y2_09192019.xlsx', sheet='Full List+Map')
-zips <- zctas(cb=TRUE, starts_with = '98')
-zips <- zips[zips$ZCTA5CE10 %in% as.character(dat$`Zip Codes`),]
+zips <- zipsAll[zipsAll$ZCTA5CE10 %in% as.character(dat$`Zip Codes`),]
 zips <- st_transform(st_as_sf(zips),4326)
 
 centers <- st_coordinates(st_centroid(zips))
@@ -34,22 +45,11 @@ zips <- cbind(zips,centers)
 geojson_write(zips, geometry = "polygon", file = 'sfs_domain_zipcodes.geojson')
 
 
-## merge neighborhoods in city with pumas in surround
-# keep state, name, lowest geoid, puma, tract, domain
-domain <- st_sf( regional_name = c(paste('Seattle--',as.character(city$NEIGHBO),sep=''), as.character(surround$NAME10)),
-                      STATE = 53,
-                      GEOID = c(rep('NA',nrow(city)), paste('53',surround$residence_puma,sep='')),
-                      domain = 'SFS_year2',
-                      geometry = c(city$geometry,surround$geometry)
-                     )
-
-
 leaflet() %>% 
   addTiles() %>% 
   addPolygons(data = domain) %>%
   addPolygons(data = st_union(zips), fillOpacity = 0, color = "red") 
 
-
 # intersect zips with domain
 # http://rpubs.com/sogletr/sf-ops
 x1 <- st_intersects(domain, st_union(zips))
@@ -62,33 +62,47 @@ x2 <- map_lgl(x1, function(x) {
 })
 ex1 <- domain[x2,]
 
-# filter 1 from tacoma that barely touchs
+# filter 1 from tacoma that barely touchs and 
+# filter 1 from Snohomish that overlaps a zip with no residents (it only covers an airport)
 ex1$regional_name
-ex2 <- ex1[-c(21),]
+ex2 <- ex1[-c(21,15),]
+ex2$regional_name <- as.character(ex2$regional_name)
+
+
+# clip rural pumas
+ex2[15,] <- st_intersection(ex2[15,],st_union(zips))
+# ex2[15,]$regional_name <- "King County (Southeast)--Maple Valley & Covington"
+
+# 6 needs to filtered to major component only
+tmp <- st_intersection(ex2[16,],st_union(zips))
+geoms <- lapply( tmp$geometry, `[` )
+tmp2<-as.data.frame(ex2[16,])
+st_geometry(tmp2) <-st_sfc(geoms[[1]][[14]])
+ex2[16,] <- tmp2
+# ex2[16,]$regional_name <- "King County (Northeast)--Cottage Lake, Union Hill & Novelty Hill"
 
 # centroids
 centers <- st_coordinates(st_centroid(ex2))
 
 # adjust centroids of 3 that are mostly rural to be closer to population-weighted
-centers[15,1] <- -122.1
-centers[16,] <- c(-122.08,47.38)
-centers[17,] <- c(-122.04,47.71)
+ex2$regional_name
+centers[14,] <- c(-122.045,47.540)
+centers[15,] <- c(-122.06,47.36)
+centers[32,] <- c(-122.25,47.3)
 
 # adjust centroid of 1 including vashon to be closer to population-weighted
-centers[18,] <- c(-122.33,47.31)
+centers[17,] <- c(-122.33,47.31)
 
 colnames(centers) <- c('lon','lat')
 ex3 <- cbind(ex2,centers)
 
-
 leaflet() %>% 
   addTiles() %>% 
   # addPolygons(data = domain) %>%
-  addPolygons(data = st_union(zips), fillOpacity = 0, color = "red") %>%
+  # addPolygons(data = st_union(zips), fillOpacity = 0, color = "red") %>%
   addPolygons(data = ex2, fillOpacity = 0, color = "green") %>%
-  addMarkers( data = ex3, lng=ex3$lon, lat=ex3$lat)
+  addMarkers( data = ex3, lng=ex3$lon, lat=ex3$lat) 
 
-ex3$regional_name <- droplevels(ex3$regional_name)
 
 geojson_write(ex3, geometry = "polygon", file = 'sfs_domain_neighborhood+puma.geojson')