Unable to intersect polygons in R even though CRS matches?

Question

I am trying to use areal interpolation to distribute the values from a source polygons layer into smaller, differently-shaped polygons from a target layer, and have been unable to do so. I have been using the sf package in in R, and specifically the aw_interpolate() function within the areal package. After encountering issues where the values of interest just were not showing up in the target layer after running the areal interpolation, all appearing as "NA", I noticed that I was not even able to intersect the polygons to begin with, leading me to speculate that these polygons layer are not even overlapping to begin with as they should be. I will explain further below.

I have the following grid of the US with varying-sized cells, each containing emissions values, which is my source layer:

I want to then take the values from this source layer and distribute them amongst the block group polygons in this target layer (yes, taking the data layer of the whole US and "downscaling" the data to the polygons for a specific county):

Of course, this means that only the grid cells that actually overlap with my county block group polygons are needed, and the rest of the grid cells in this state, and the rest of the US, are ignored.

However, running the areal interpolation yields all NA values in my target block group polygons for my field of interest. I wanted to break apart the process, as I was suspicious about the layers somehow not overlapping, and tried to run a simple intersection in R, just to make sure this actually worked.

Here is my R code:

install.packages("sf")
install.packages("areal")
library(sf)
library(areal)

US_emissions <- st_read("C:/Users/Name/Desktop/Map/US_Output/US_Output_.shp")

County_blockgroups <- st_read("C:/Users/Name/Desktop/Map/County_Blockgroups/County_Blockgroups.shp")

#Get these two layers in the same CRS
st_crs(US_emissions) <- st_crs(County_blockgroups)
US_emissions <- st_transform(US_emissions, crs = st_crs(County_blockgroups))

And then just to make sure my two layers are in the same CRS, I then run:

US_emissions

and receive:

Simple feature collection with 46998 features and 20 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -2736000 ymin: -2088000 xmax: 2448000 ymax: 1944000
Projected CRS: NAD83 / California zone 4 (ftUS)
First 10 features:

and then run:

County_blockgroups

and receive:

Simple feature collection with 589 features and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 5999505 ymin: 1851910 xmax: 6749006 ymax: 2460485
Projected CRS: NAD83 / California zone 4 (ftUS)
First 10 features:

Great, it looks like these two layers are in the same CRS and should overlap fine!

So then I try to intersect the two layers, running the following:

intersect <- st_intersection(County_blockgroups, US_emissions)
intersect

and receive:

Warning: attribute variables are assumed to be spatially constant throughout all geometries
Simple feature collection with 0 features and 36 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: NAD83 / California zone 4 (ftUS)
 [1] AWATER10   TRACTCE10  STLength   STArea     ALAND10    BLKGRPCE10 FUNCSTAT10
 [8] ID         NAMELSAD10 COF_Census COUNTYFP10 STATEFP10  MTFCC10    GEOID10   
[15] INTPTLAT10 INTPTLON10 AsianD     BasePM25   BlackD     LatinoD    NH3       
[22] NOx        NativeD    PNH4       PNO3       PSO4       PrimPM25   SOA       
[29] SOx        TotalPM25  TotalPop   TotalPopD  WhitNoLatD WindSpeed  PolygonID 
[36] Area_Sourc geometry  
<0 rows> (or 0-length row.names)

It looks like the intersection "worked", but then again, it looks like nothing actually intersected, seeing that there are 0 features and 0 rows produced. If the intersection does not work, then I would think that this is the reason the areal interpolation did not work, since to my understanding, intersection is a needed step in areal interpolation. Though I do not understand why the intersection did not work if my source layer and target layer are both successfully read into R, and both appear to be in the same CRS. Where might I be going wrong here?

Answer 1

This bit is possibly wrong:

#Get these two layers in the same CRS
st_crs(US_emissions) <- st_crs(County_blockgroups)
US_emissions <- st_transform(US_emissions, crs = st_crs(County_blockgroups))

When you first create US_emissions using st_read it should get a CRS from the shapefile. It should be the correct CRS of the data. That first line above - st_crs(US_emissions) <- st_crs(County_blockgroups) will overwrite the CRS of US_emissions with - possibly - the wrong CRS!

Only the second of those lines should be necessary - it computes the transformed coordinate numbers of US_emissions into the CRS of County_blockgroups and returnes an object with that CRS code assigned to it.

You should always check that the bounding box and the st_crs of spatial objects are correct when you read them, and you should also show this for data read into questions here (especially if you cant supply the data).