Search code examples
rspatialshapefile

subsetting zipcodes from shapefile and plotting in R

I have a Census Shapefile of the United States by ZCTA. Am only interested in plotting NY. I can read in the whole shapefile and plot only NY by limiting the longitude and latitude but it takes a long time and R keeps crashing.

Instead, I'd like to reduce to only zipcodes in NY. I'm reading in the shapefile, converting to dataframe and then subsetting to the zipcodes/ZCTAs in NY. When I plot the result, however, it's not a full map. Is this the right approach? My R code and plot are below:

### read in the ZCTA shapefile for US
us_modzcta <- "Data/us_zcta/tl_2018_us_zcta510.shp"
us <- shapefile(us_modzcta)

## convert to dataframe 
us_df <- fortify(us,region="ZCTA5CE10")

## subset - keep only zip codes in NY - starts with 10, 11, 12, 13, 14
ny_df <- us_df[which(startsWith(us_df$id,c("11","12","13","10","14"))),]
p<- plot(ny_df)

enter image description here

Solution

  • Why so complicated path via raster and ggplot2? Please have a look on sf package:

    library(sf)
library(tidyverse)
us_modzcta <- "census/tl_2018_us_zcta510.shp"
ny <- st_read(us_modzcta) |>
  filter(grepl("^1[0-4]", ZCTA5CE10))

plot(ny$geometry)

    Created on 2022-03-05 by the reprex package (v2.0.1)