Baltimore.md

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)

#   rw <- data.frame( ## Simplified neighborhood rectangle
#      Longitude = c(-8528150, -8528500, -8528500, -8528150),
#      Latitude  = c( 4771475,  4771475,  4771880,  4771880)
# )
  
prj <- "EPSG:3857"
lat_max <- 39.3525
long_max <- -76.617
lat_min <- 39.3455
long_min <- -76.6095
rw2 <- terra::project(cbind(c(long_min, long_max), c(lat_min, lat_max)), to = prj, from = "OGC:CRS84")
  
## perhaps a 7 got in there in place of the 4 in the y coord?
# rw2
#         [,1]    [,2]
#[1,] -8528131 4744189
#[2,] -8528965 4745192


library(vapour)
library(ximage)  ## hypertidy/ximage
im <- gdal_raster_image(dsn::wms_virtualearth_street(), target_ext = c(range(rw2[,1]), range(rw2[,2])) + c(-1, 1, -1, 1) * 1e2, 
                        target_crs = prj, target_dim = dev.size("px"), resample = "cubic")
#> Warning in warp_general_cpp(dsn, target_crs, target_ext, target_dim,
#> target_res, : no source crs, target crs is ignored

op <- par(mar = rep(0, 4))
ximage(im, asp = 1)
abline(v = range(rw2[,1]), h = range(rw2[,2]))

^{Created on 2023-06-13 by the reprex package (v2.0.1)}

and apologies, I don't know how 3395 got in there I think 3857 is the right code to use given the the tile server (not sure it would matter substantially though).

Mike, thanks, again, for all your help. You got me started down the correct path.

Here's what's finally working for me:

library(tidyverse)
library(tidycensus)
library(sf)
library(tmap)
library(tigris)
options(tigris_use_cache = TRUE)
library(tmaptools)

rw_block_list <- c("3000", "3001", "3002", "3005", "3006", "3007",
                   "3008", "3009", "3010", "3011", "3012")

## Get the RW blocks from the census:
rw_blocks <- blocks(state = "MD",
                    county = "Baltimore city",
                    year = "2020") %>%
    filter(substr(GEOID20, 6, 11) == "271101" &
           substr(GEOID20, 12, 15) %in% rw_block_list)

## Create a map of just the RW blocks:
rw_base_blocks <- read_osm(bb(rw_blocks, ext = 1.3))

tmap_mode("plot")

## Line below gives map in meters
(RW_block_map <- tm_shape(rw_base_blocks, projection = 6487) +
## Line below gives map in degrees
## (RW_block_map <- tm_shape(rw_base_blocks, projection = 6487) +
     tm_rgb() +
     tm_shape(rw_blocks) +
     tm_fill("MAP_COLORS", alpha = 0.2, palette = "Accent") +
     tm_borders() +
     tm_scale_bar() +
     ## tm_grid() + tm_xlab("Long") + tm_ylab("Lat") +
     tm_grid() +
     tm_layout(title = "Radnor-Winston Neighborhood")
)

tmap_save(RW_block_map, "rw_map.png")

=======================================

Based on the map in meters, I can pick off the coordinates of a polygon in my neighborhood, and convert it to the correct location in degrees:

===========================================

## Test polygon:
base_x <- 433000
base_y <- 186000
rw_neigh_pg_m <- data.frame(
    matrix(
        c(300, 900,
          500, 900,
          500, 700,
          300, 700
        ),
        ncol = 2, byrow = TRUE)
) %>% + matrix(c(rep(base_x, nrow(.)), rep(base_y, nrow(.))),
               nrow = nrow(.)) %>%
sf::st_as_sf(coords = c(1,2), dim = "XY") %>%
summarize(geometry = st_combine(geometry)) %>%
st_cast("POLYGON") %>%
st_set_crs(6487)

## Convert to degrees:
(rw_neigh_pg_d <- rw_neigh_pg_m %>%
    st_transform(4326)
)

=============================================

Thanks, again, for all your help. I would have spent many additional hours on randomly trying things without your guidance.

-Kevin