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

Kevin Zembower kev|n @end|ng |rom zembower@org
Tue Jun 13 16:40:21 CEST 2023


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 using zembower.org 
> <mailto:kevin using 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
>     <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>
>     <mailto: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>
>     <mailto: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>
>      >     <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>
>     <mailto:mdsumner using gmail.com <mailto:mdsumner using gmail.com>>
> 
> 
> 
> -- 
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsumner using gmail.com <mailto:mdsumner using gmail.com>




More information about the R-sig-Geo mailing list