[R-sig-Geo] Mapping my own polygon?

[Kevin Zembower]

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 using r-project.org <mailto:r-sig-geo using 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
> 
>     _______________________________________________
>     R-sig-Geo mailing list
>     R-sig-Geo using r-project.org <mailto:R-sig-Geo using r-project.org>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> 
> 
> 
> -- 
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsumner using gmail.com <mailto:mdsumner using gmail.com>