[R-sig-Geo] Mapping my own polygon?
Michael Sumner
md@umner @end|ng |rom gm@||@com
Tue Jun 13 03:47:07 CEST 2023
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 <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/, 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>
>
>
--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner using gmail.com
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list