Mapping my own polygon?

Hi, Mike, thanks for answering me from across the world!

I created my polygon of my neighborhood the old fashion way: by printing
a map with a grid in meters, and using a pair of navigational dividers
to pick off the x and y coordinates of the borders of my neighborhood.
You can see a map of my neighborhood, drawn by someone else, at
https://www.radnorwinston.org/, if you're curious.

Here's code to generate a map of this area, one with degrees lat/long,
and the one I used, with a grid of meters:

===========================================
## This gives degrees lat long:
library(tidyverse)
library(tigris)
options(tigris_use_cache = TRUE)
library(sf)
library(OpenStreetMap)

lat_max <- 39.3525
long_max <- -76.617
lat_min <- 39.3455
long_min <- -76.6095
nw <- c(lat_max, long_max)
se <- c(lat_min, long_min)

rw_map <- openmap(nw, se,
                   type = "osm",
                   mergeTiles = TRUE) %>%
     openproj() %>%
     OpenStreetMap::autoplot.OpenStreetMap() +
     xlab("long") + ylab("lat")

rw_map

## This gives map with grid in meters:
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")

(RW_block_map <- tm_shape(rw_base_blocks) +
      tm_rgb() +
      tm_shape(rw_blocks) +
      tm_fill("MAP_COLORS", alpha = 0.2, palette = "Accent", n = 10) +
      tm_borders() +
      tm_scale_bar() +
      tm_grid() + tm_xlab("Long") + tm_ylab("Lat") +
      tm_layout(title = "Radnor-Winston Neighborhood")
)
===================================================

You're right, that I did pick my example "1 POINT (-76.61246 39.35010)"
just by picking a random point in my neighborhood, using Google Maps.

In the second example that gives the grid in meters, I just assumed that
this was CRS 6487, due to it being Maryland and in meters. Was this a
bad assumption? How can I tell what the CRS of the return of the call to
"read_osm(bb(rw_blocks, ext = 1.3))" is? When I display it, it says:

  > rw_base_blocks

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
    Min. 1st Qu. Median    Mean 3rd Qu. Max.
X     0     217    224 219.915     238  255
dimension(s):
      from  to   offset    delta                   refsys
values x/y
x       1 663 -8528854  1.19942 WGS 84 / Pseudo-Mercator
NULL [x]
y       1 907  4772338 -1.19981 WGS 84 / Pseudo-Mercator
NULL [y]
band    1   3       NA       NA                       NA red  , green,
blue

 >

Thanks, again, for your work to help me. Thanks also for introducing me
to Github Gist, which I had never heard of before. I'm going to go to
that page next and see if I can put this response in there.

Take care.

-Kevin

On 6/11/23 04:04, Michael Sumner wrote:

I took a guess that your coordinates are not in EPSG:6487 but in global
Mercator (EPSG:3857) which seems to give a reasonable region from online
image servers.

https://gist.github.com/mdsumner/e6997c2f4a54c743e078aca8401537a0
<https://gist.github.com/mdsumner/e6997c2f4a54c743e078aca8401537a0>

If that looks ok?     Then, simply replace 6487 with 3857 in your code,
but also please take care to track down where your neighbourhood
coordinates for the x,y range came from.  3857 is the infamous global
google mercator, and that might be what you're using or simply close to
it, you should make sure you know. :)

I'm using in-dev code in the example, it's not because I want you to
use  that or advise you to - it's just to keep a record of what I did. I
experimented with the code and scale to guess at what might be the

problem.

If the pic in my gist is not on the right track I'm happy to follow up,
best of luck!

HTH, Mike



On Sun, Jun 11, 2023 at 5:27?AM Kevin Zembower via R-sig-Geo
<r-sig-geo at r-project.org <mailto:r-sig-geo at r-project.org>> wrote:

    In my continuing work on reporting on US Census data for my
    neighborhood, I'd like to draw a map of the boundaries of it. I was
    successful in creating and printing an OSM basemap, with the US

Census

    blocks that make up my neighborhood on it.

    Now, I'd like to create my own polygon, of the boundaries of my
    neighborhood, because the census blocks don't line up exactly with

the

    neighborhood boundaries. I need help creating a polygon that I can
    submit to read_osm() that will correctly return an OSM map of my

area.

    Here's what I've tried so far:
    ## Reproducible simple example:
    library(tidyverse)
    library(sf)

    rw <- data.frame( ## Simplified neighborhood rectangle
          Longitude = c(-8528150, -8528500, -8528500, -8528150),
          Latitude  = c( 4771475,  4771475,  4771880,  4771880)
    )

    rw ## Returns (as expected):

      > rw

        Longitude Latitude
    1  -8528150  4771475
    2  -8528500  4771475
    3  -8528500  4771880
    4  -8528150  4771880

      >

    rw %>%
          st_as_sf(coords = c("Longitude", "Latitude"), dim = "XY") %>%
          st_set_crs(6487) %>% ## CRS 6487 is NAD83 (2011) Maryland in
    meters
          st_transform(crs = 4269) ## CRS 4269 is NAD83

    ## Returns:

    Simple feature collection with 4 features and 0 fields
    Geometry type: POINT
    Dimension:     XY
    Bounding box:  xmin: 171.777 ymin: 24.65904 xmax: 171.7818 ymax:
    24.66314
    Geodetic CRS:  NAD83
                         geometry
    1 POINT (171.7818 24.66192)
    2 POINT (171.7806 24.65904)
    3  POINT (171.777 24.66026)
    4 POINT (171.7781 24.66314)

      >

    I expected the POINTS to look like:
    1 POINT (-76.61246 39.35010)

    Can anyone suggest what I'm doing wrong? Thanks so much in advance.
    I've
    worked on it all day today, without making much progress.

    -Kevin