Do a point in polygon join using sf package to get polygon attributes where points land - in this case the NHDPlus catchment ID where points land

SF_Point_Poly.R

library(sf)
gages <- st_read('C:/Users/mweber/Temp/gages_test.shp')
class(gages)

# Let's pretend it's a csv file we read in rather than a shapefile - I'll strip out the 'spatial' 
# part (the geometry column) and make it just a data.frame with Lon and Lat columns, then show
# how to promote it back

st_geometry(gages) <- NULL
class(gages)
# now imagine we've read in a flat file with coordinate info to a data frame 'gages' - since we're just using the 'Lat'
# and 'Lon' columns, we apply NAD83 as the coordinate ref system rather than projection the shapefile was in
gages <- st_as_sf(gages, coords = c("LON_SITE", "LAT_SITE"), crs = 4269,agr = "constant")
class(gages)

# These are gages in Oregon so I'm reading in NHDPlus catchments for hydro-region 17
NHD_cats <- st_read('H:/NHDPlusV21/NHDPlusPN/NHDPlus17/NHDPlusCatchment/Catchment.shp')
# Projection have to be the same to join features - check
st_crs(NHD_cats) == st_crs(gages)
# Just encoded slightly differently
gages <- st_transform(gages, st_crs(NHD_cats))

# Now do spatial join of points and polys
z <- st_join(gages, NHD_cats)
head(z)