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 <ke...@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@r-project.org <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: > > ## 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 <mailto:R-sig-Geo@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: mdsum...@gmail.com <mailto:mdsum...@gmail.com> > > -- 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
