Re: [R-sig-Geo] 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
>
>
>
>
>
>
>
> On Mon, Jun 12, 2023 at 12:58 AM Kevin Zembower > 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/
Re: [R-sig-Geo] Mapping my own polygon?
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
On Mon, Jun 12, 2023 at 12:58 AM Kevin Zembower 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/, 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. MedianMean 3rd Qu. Max.
> X 0 217224 219.915 238 255
> dimension(s):
> from to offsetdelta 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]
> band1 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
> >
> >
> > 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
> > mailto:r-sig-geo@r-project.or
Re: [R-sig-Geo] 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. MedianMean 3rd Qu. Max.
X 0 217224 219.915 238 255
dimension(s):
from to offsetdelta 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]
band1 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
>
>
> 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
> mailto:r-sig-geo@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:
> ## Reproducib
Re: [R-sig-Geo] Mapping my own polygon?
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
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@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@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@gmail.com
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
