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
