0

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

Michelle
  • 55
  • 1
  • 7
  • Hi Michelle, it would be easier to help you if you could make your example reproducible. For this, it would be helpful to make accessible your shapefile `tl_2018_us_zcta510.shp`. Cheers. – lovalery Mar 04 '22 at 20:46
  • 1
    can be downloaded here: https://www.census.gov/cgi-bin/geo/shapefiles/index.php . select 2018 as year and Zip Code tabulation areas as "Layer Type" Also tried subsetting by latitude and longitude and that produces a similar plot – Michelle Mar 04 '22 at 20:59
  • I'm really sorry: I manage to download the data but I can't unzip them (my pc crashes). Therefore it is difficult for me to help you. Maybe someone else will not have the same issues and will be able to provide the help you are looking for. – lovalery Mar 04 '22 at 21:45
  • See https://stackoverflow.com/questions/13443372/simple-way-to-subset-spatialpolygonsdataframe-i-e-delete-polygons-by-attribut perhaps ny_df <- us_df@data[us_df@data$id %in% c(("11","12","13","10","14")] – Susan Switzer Mar 04 '22 at 22:53
  • @SusanSwitzer tried your suggestion but couldn't get it to work – Michelle Mar 06 '22 at 19:38

1 Answers1

1

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)

Grzegorz Sapijaszko
  • 1,913
  • 1
  • 5
  • 12