0

I have two datasets.

NY Zipcodes from that website: https://data.cityofnewyork.us/Business/Zip-Code-Boundaries/i8iw-xf4u/data?no_mobile=true

zc <-  st_read("ZIP_CODE_040114/ZIP_CODE_040114.shp")

And the other is the white population of New York County with get_acs

library(tidycensus)
library(tidyverse)
census_api_key("mykey")
ny_wp <- get_acs(geography = "tract", variables = "B03002_003",
                state = "NY", county = "New York", geometry = TRUE, year = 2020)

now I did these steps to do an st_join

zc <- st_set_crs(zc, 2263)
zc <- st_transform(zc,crs= 2263) 


ny_wp  <- st_set_crs(ny_wp , 2263)
ny_wp  <- st_transform(ny_wp ,crs= 4326) 
ny_wp  <- st_transform(ny_wp , crs = st_crs(zc))

st_join(zc,ny_wp)

But I have no numbers for white population. I want the white population per zipcode?

1 Answer 1

0

Why do this?

zc <- st_set_crs(zc, 2263)
zc <- st_transform(zc,crs= 2263) 

from the start, zc already is epsg 2263:

> zc = st_read("./ZIP_CODE_040114.shp", quiet=TRUE)
> st_crs(zc)$epsg
[1] 2263

so the first line is telling sf what it already knows, and the second line is transforming it to the coordinate system that it already is.

Then why do this?

ny_wp  <- st_set_crs(ny_wp , 2263)
ny_wp  <- st_transform(ny_wp ,crs= 4326) 
ny_wp  <- st_transform(ny_wp , crs = st_crs(zc))

ny_wp comes as epsg:4269

> ny_wp <- get_acs(geography = "tract", variables = "B03002_003",state = "NY", county = "New York", geometry = TRUE, year = 2020) 
Getting data from the 2016-2020 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
Using FIPS code '36' for state 'NY'
Using FIPS code '061' for 'New York County'
> st_crs(ny_wp)$epsg
[1] 4269

...so your first line is setting it to the wrong value. This breaks all the subsequent transforms because you've told sf its in epsg:2263 and it isnt, its epsg:4269. I'm not sure then why you do transforms when you can do that in one, but no matter, its already broken by setting the CRS to the wrong value.

  • If a geographic data set has a correct CRS, then don't use st_set_crs.
  • To transform one data set to the coords of another, use one st_transform.
     > zc = st_read("./ZIP_CODE_040114.shp", quiet=TRUE)
     > ny_wp <- get_acs(geography = "tract", variables = "B03002_003",state = "NY", county = "New York", geometry = TRUE, year = 2020) 
     Getting data from the 2016-2020 5-year ACS
     Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
     Using FIPS code '36' for state 'NY'
     Using FIPS code '061' for 'New York County'
     > zc = st_transform(zc, st_crs(ny_wp))
     > 
     > jj = st_join(zc, ny_wp)

Assuming estimate is the value of the variable, and B03002_003 is the code for the variable of interest, you've got it:

> jj[200:210,c("ZIPCODE","variable","estimate")]
Simple feature collection with 11 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -73.96108 ymin: 40.79663 xmax: -73.90219 ymax: 40.81427
Geodetic CRS:  NAD83
First 10 features:
     ZIPCODE   variable estimate                       geometry
48.1   10454 B03002_003      290 POLYGON ((-73.93176 40.8079...
48.2   10454 B03002_003      308 POLYGON ((-73.93176 40.8079...
49     10026 B03002_003     1034 POLYGON ((-73.9446 40.80323...
49.1   10026 B03002_003     1433 POLYGON ((-73.9446 40.80323...
49.2   10026 B03002_003     2858 POLYGON ((-73.9446 40.80323...

which you can plot:

plot(jj[,"estimate"])

enter image description here

I notice this is overlaying zip code areas with census tracts so more careful consideration is needed to the messy intersection of these two differing polygon sets. But I've hopefully explained your coordinate system issues.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.