1

I am trying to perform a spatial join of two sf shape files.I am losing all information from the second data set (i.e output_inmap). Whichever dataset is placed second will return all NA values. Anyone know what could be happening?

output_inmap <- st_read("processed/ceidars_data_inmap.shp")
output_inmap <-st_transform(output_inmap, crs=3310)

unzip("census-tract.zip")
census_tracts <- st_read("census-tract/tl_2019_06_tract.shp")

st_transform(census_tracts, crs = 3310)
st_transform(output_inmap, crs = 3310)

TC_1<- st_join(census_tracts, output_inmap)

I am losing all information from the second data set (i.e output_inmap). Whichever dataset is placed second will return all NA values. Anyone know what could be happening?

Tyler Cobian
  • 21
  • 2
  • Hard to say unless you post a reproducible example: https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-examplehttps://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – Bill O'Brien Aug 05 '20 at 19:10
  • What are the types of the datasets? st_join works if both classes are sf objects. I use geo_join {tigris} to join spatial data with dataframes. – Susan Switzer Aug 05 '20 at 20:34
  • @SusanSwitzer, they are both sf objects, but I will try geo_join() – Tyler Cobian Aug 05 '20 at 21:33
  • 1
    @BillO'Brien, you are right sorry, I am transitioning from being a passive user of Stack to an active member. The question is updated with more information. – Tyler Cobian Aug 05 '20 at 21:45
  • Not sure what's happening there, but `st_transform()` does not modify anything in place, so you need to write back to the original variable. Also maybe your shape files are in fact non-overlapping. Can you check for validity of your shapes `sf::st_is_valid()`? – dmi3kno Aug 05 '20 at 21:54
  • @dmi3kno, st_is_valid() is coming up as valid. I tried to specify the crs in the st_read() but still having the same issue ```output_inmap <- st_read("processed/ceidars_data_inmap.shp", crs = 4269) census_tracts <- st_read("census-tract/tl_2019_06_tract.shp", crs = 4269)``` – Tyler Cobian Aug 05 '20 at 22:15

1 Answers1

0

Your second st_transform (of the census tracts) seems to be leading nowhere; consider this code (slightly adjusted via dplyr style pipe) to ensure both spatial objects are on the same CRS.

You may also consider setting parameter left of the sf::st_join() call (by default true) to false = change behaviour from left (preserving) to inner (filtering) style join. Sometimes this makes for a more concise code.

library(sf)
library(dplyr)

output_inmap <- st_read("processed/ceidars_data_inmap.shp") %>%
   st_transform(crs=3310)

unzip("census-tract.zip")
census_tracts <- st_read("census-tract/tl_2019_06_tract.shp") %>%
   st_transform(crs = 3310)


TC_1<- st_join(census_tracts, output_inmap)
Jindra Lacko
  • 7,814
  • 3
  • 22
  • 44
  • I tried that code and I am still loosing data from the second data input in the st_join. – Tyler Cobian Aug 12 '20 at 15:48
  • `STATEFP COUNTYFP TRACTCE GEOID NAME NAMELSAD MTFCC 1 06 037 139301 06037139301 1393.01 Census Tract 1393.01 G5020 FUNCSTAT ALAND AWATER INTPTLAT INTPTLON field_1 CARB_ID latitude 1 S 2865657 0 +34.1781538 -118.5581265 NA NA NA longitude NH3 PM2_5 NOx SOx VOC height diam temp velocity 1 NA NA NA NA NA NA NA NA NA NA geometry 1 MULTIPOLYGON (((-118.5715 3...` – Tyler Cobian Aug 12 '20 at 15:57
  • And are you 100% certain your objects are overlapping? What does `mapview::mapview(census_tracts) + mapview::mapview(output_inmap)` show? – Jindra Lacko Aug 13 '20 at 20:24