What's the Zip Code?

library(tmaptools)
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tmap)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

2 Addresses without zip codes

address<-c("630 Old Country Rd, Garden City, NY ", "20 W 34th St., New York, NY")

#OSM=Open Street Map
geocoded_addresses<-geocode_OSM(address)
print(geocoded_addresses)
##                                  query      lat       lon  lat_min  lat_max
## 1 630 Old Country Rd, Garden City, NY  40.73967 -73.61354 40.73962 40.73972
## 2          20 W 34th St., New York, NY 40.74865 -73.98530 40.74860 40.74870
##     lon_min   lon_max
## 1 -73.61359 -73.61349
## 2 -73.98535 -73.98525

Dataset is from censusreporter.org

#Zip file from censusreporter.org of NY State Total Population as a shape file
unzip("/Users/Sangeetha/Downloads/acs2020_5yr_B01003_86000US14742.zip", exdir = "ny_zip_shapefile", junkpaths = TRUE, overwrite = TRUE)

Reading in data & plotting a map

zipcode_geo <- st_read("/Users/Sangeetha/ny_zip_shapefile/acs2020_5yr_B01003_86000US14742.shp")
## Reading layer `acs2020_5yr_B01003_86000US14742' from data source 
##   `/Users/Sangeetha/ny_zip_shapefile/acs2020_5yr_B01003_86000US14742.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1786 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -79.76259 ymin: 40.47658 xmax: -71.77749 ymax: 45.01586
## Geodetic CRS:  WGS 84
zipcode_geo <-filter(zipcode_geo, name!="New York")
#quick theme map from tmaptools
qtm(zipcode_geo, title="NY Population Chloropleth Map") +  tm_legend(show = FALSE) 

#crs =coordinate reference systems
#dataframe converted into sf geospatial object
point_geo <-st_as_sf(geocoded_addresses, coords=c(x="lon", y="lat"), crs=4326)

my_results<-st_join(point_geo, zipcode_geo, join=st_within)
#st_join(point_sf_object, polygon_sf_object, join=join_type)
print(my_results)
## Simple feature collection with 2 features and 9 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.9853 ymin: 40.73967 xmax: -73.61354 ymax: 40.74865
## Geodetic CRS:  WGS 84
##                                  query  lat_min  lat_max   lon_min   lon_max
## 1 630 Old Country Rd, Garden City, NY  40.73962 40.73972 -73.61359 -73.61349
## 2          20 W 34th St., New York, NY 40.74860 40.74870 -73.98535 -73.98525
##          geoid  name B01003001 B01003001e                   geometry
## 1 86000US11530 11530     27669        543 POINT (-73.61354 40.73967)
## 2 86000US10001 10001     25026       1759  POINT (-73.9853 40.74865)

Reprinting with my selected columns

print(my_results[, c("query", "name","geometry")])
## Simple feature collection with 2 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.9853 ymin: 40.73967 xmax: -73.61354 ymax: 40.74865
## Geodetic CRS:  WGS 84
##                                  query  name                   geometry
## 1 630 Old Country Rd, Garden City, NY  11530 POINT (-73.61354 40.73967)
## 2          20 W 34th St., New York, NY 10001  POINT (-73.9853 40.74865)

Plotting my 2 points on a map

tm_shape(zipcode_geo)+
  tm_fill()+
  tm_shape(my_results)+
  tm_bubbles(col="blue", size =0.25)