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> _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo