Mapping my own polygon?
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
On 6/12/23 21:47, Michael Sumner wrote:
The OpenStreetMap package is using PseudoMercator (quite common for image tile servers) which is EPSG:3857. Using your long/lat min max the numbers give the right neighbourhood, perhaps a transcribe error in the original 'rw' data frame numbers for Y? https://gist.github.com/mdsumner/e6997c2f4a54c743e078aca8401537a0?permalink_comment_id=4597644#gistcomment-4597644 <https://gist.github.com/mdsumner/e6997c2f4a54c743e078aca8401537a0?permalink_comment_id=4597644#gistcomment-4597644> On Mon, Jun 12, 2023 at 12:58?AM Kevin Zembower <kevin at zembower.org <mailto:kevin at zembower.org>> wrote: 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/ <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 <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>
<mailto: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
>
>? ? ?_______________________________________________
>? ? ?R-sig-Geo mailing list
> R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>
<mailto:R-sig-Geo at r-project.org <mailto:R-sig-Geo at r-project.org>>
>
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsumner at gmail.com <mailto:mdsumner at gmail.com>
<mailto:mdsumner at gmail.com <mailto:mdsumner at gmail.com>> -- Michael Sumner Software and Database Engineer Australian Antarctic Division Hobart, Australia e-mail: mdsumner at gmail.com <mailto:mdsumner at gmail.com>