https://stat.ethz.ch/pipermail/r-sig-geo/2023-June/029294.html
rw <- data.frame( ## Simplified neighborhood rectangle
Longitude = c(-8528150, -8528500, -8528500, -8528150),
Latitude = c( 4771475, 4771475, 4771880, 4771880)
)
prj <- "EPSG:3857"
library(vapour)
library(ximage) ## hypertidy/ximage
library(dsn) ## hypertidy/dsn
im <- gdal_raster_image(dsn::wms_virtualearth_street(), target_ext = c(range(rw$Longitude), range(rw$Latitude)) + c(-1, 1, -1, 1) * 5e4,
target_crs = prj, target_dim = dev.size("px"), resample = "cubic")
op <- par(mar = rep(0, 4))
ximage(im, asp = 1)
abline(v = rw$Longitude, h = rw$Latitude)
par(op)
library(sf)
rw %>%
st_as_sf(coords = c("Longitude", "Latitude"), dim = "XY") %>%
st_set_crs(3857) %>%
st_transform(crs = 4269) ## CRS 4269 is NAD83
Simple feature collection with 4 features and 0 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -76.61282 ymin: 39.34683 xmax: -76.60968 ymax: 39.34965
Geodetic CRS: NAD83
geometry
1 POINT (-76.60968 39.34683)
2 POINT (-76.61282 39.34683)
3 POINT (-76.61282 39.34965)
4 POINT (-76.60968 39.34965)
fwiw, I feel like it's still not the right interpretation because the actual region doesn't seem to correspond to anything in particular:
par(op)
Created on 2023-06-11 by the reprex package (v2.0.1)