Fast Reverse Geocoding Using R and a Local PostGIS Database

library(tidyverse)
library(RPostgreSQL)
library(sf)

con <- DBI::dbConnect(RPostgres::Postgres(), dbname='gis', host='localhost', 
                      port=5432, user='postgres', password='xxxxx')

# function that uses PostGIS query to return the closest
# point in the SAM database
reverseGeocode <- function(geom) {
  # convert sf geom object to text
  geom <- st_as_text(geom)
  
  # query
  sql <- str_c("
SELECT 
  *,
  \"Shape\" <-> 'SRID=26918;", geom, "'::geometry AS dist
FROM
  nys_sam_points
ORDER BY
  dist
LIMIT 1;")
  
  query <- dbGetQuery(con, sql)
  
  str_c(query$addresslabel, 
        ', ',
        query$citytownname,
        ', NY ',
        query$zipcode
        ) %>% return
}

# obtain centroids for each county
ny.county.centroids <- tigris::county_subdivisions('ny') %>% 
  st_centroid()

# get addresses of properties by reverse gecoding
ny.county.centroids %>%
  st_transform(26918) %>%   # use same projection 
                            # as your postgis database
  mutate(
    # obtain address using purr and map_chr using the 
    # reverseGeocode Function defined above
    Address = map_chr(geometry, reverseGeocode)
  )

library(tidyverse)
library(RPostgreSQL)
library(sf)

con <- DBI::dbConnect(RPostgres::Postgres(), dbname='gis', host='localhost', 
                      port=5432, user='postgres', password='xxxxx')

# function that uses PostGIS query to interpolate an address from the 
# streeet high and low adress
reverseGeocode <- function(geom) {
  # convert sf geom object to text
  geom <- st_as_text(geom)
  
  # query
  sql <- str_c("
SELECT 
  CAST (
  (leftfromaddress +
  (righttoaddress - leftfromaddress) *  ST_LineLocatePoint(ST_LineMerge(roadline), 'SRID=26918;", geom, "'::geometry)) AS integer) AS hno,
  completestreetname, leftzipname
  
  FROM (
SELECT 
  CAST(leftfromaddress AS double precision), CAST(righttoaddress AS double precision), completestreetname, leftzipname,
  \"Shape\" as roadline,
  \"Shape\" <-> 'SRID=26918;", geom, "'::geometry AS dist
FROM
  nys_sam_streets
ORDER BY
  dist
LIMIT 1)")
  
  
  dbGetQuery(con, sql)
  
  query <- dbGetQuery(con, sql)
  
  str_c(query$hno %>% as.character() %>% replace_na(''), 
        ' ',
        query$completestreetname,
        ', ',
        query$leftzipname,
        ', NY ',
        query$leftpostal
  ) %>% return
}

# obtain centroids for each county
ny.county.centroids <- tigris::county_subdivisions('ny') %>% 
  st_centroid()

# get addresses of properties by reverse gecoding
ny.county.centroids %>%
  st_transform(26918) %>%   # use same projection 
                            # as your postgis database
  mutate(
    # obtain address using purr and map_chr using the 
    # reverseGeocode Function defined above
    Address = map_chr(geometry, reverseGeocode)
  )