2

I have a large set of coordinates from the critical and endangered habit federal registry. I'm trying to digitize these maps for analysis. Here's a sample of the data as an example.

library(tidyverse)
library(mapview)
library(sf)
library(ggplot2)

lat <- c(300786, 301348,  301467, 302384, 304644, 304747, 304863, 304882)
lon <- c(4215918, 4215650, 4215784, 4216077, 4211303, 4211336, 4211395, 4211457)
geodeticDA <- c("NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83")
utmZone <- c("11N", "11N", "11N", "11N", "11N", "11N", "11N", "11N")
ID <- c("A", "A", "A", "A", "B", "B", "B", "B")

df <- data.frame(lat, lon, geodeticDA, utmZone, ID)

I've successfully been able to turn the df into an sf and graph it using the following

plot_local1 <- st_as_sf(df, 
                        coords = c("lat", "lon"),
                         crs = "+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")

mapview(plot_local1, map.types = "Esri.NatGeoWorldMap") 

mapview coordinates plot

My main goal is to create two different polygons based on their ID (A or B) but I'm struggling.

I can turn them all into one polygon with this:

polygon <- plot_local1 %>% 
# dplyr::group_by(ID) %>% 
  dplyr::summarise() %>%
  st_cast("POLYGON")

mapview(polygon, map.types = "Esri.NatGeoWorldMap") 

Single polygon

However, if I uncomment the group_by(), as is suggested by other sites/questions, I get this error

Error: ! All columns in a tibble must be vectors. x Column geometry is a sfc_POINT/sfc object.

I've managed to create multiple polygons with the following code but when I graph it it's not correct.

polys <- st_sf(
  aggregate(
    plot_local1$geometry,
    list(plot_local1$ID),
    function(g){
       st_cast(st_combine(g),"POLYGON")
    }
   ))

mapview(polys, map.types = "Esri.NatGeoWorldMap") 

multiple polygon attempt

Certainly I could create separate data frames for each polygon and then combine them together later but given the amount of data I'm working with that really isn't practical. I've also updated all of the packages I'm using so I know that's not the issue. Any and all help is appreciated!

  • 1
    `sfheaders::sf_polygon(obj = df, x = "X", y = "Y", polygon_id = "ID")` – SymbolixAU Mar 02 '22 at 23:32
  • did my comment work for you? – SymbolixAU Mar 03 '22 at 00:51
  • @SymbolixAU I tried and I got this error Error: C stack usage 7953936 is too close to the limit. Trying to figure out how to deal with it. – JulietheFoodie Mar 03 '22 at 08:46
  • `polygon <- plot_local1 %>% group_by(ID) %>% summarize() %>% st_cast(., "POLYGON")` works for me. When running `mapview(polygon, map.types = "Esri.NatGeoWorldMap")`, it returns the last map of your post. I don't understand why this map would not be good: it seems to me consistent with the coordinates of the points. What map do you expect? Cheers. – lovalery Mar 03 '22 at 12:49
  • 1
    @lovalery When I say "wrong" I mean I can't get any actual map backgound on it. It's all just grey space. I tried your code and I got the same error as before ( All columns in a tibble must be vectors. x Column `geometry` is a `sfc_POINT/sfc` object. ) Seems like there's probably something wrong with a package or something :( – JulietheFoodie Mar 05 '22 at 01:02

1 Answers1

2

As a follow-up to your comment, I have prepared a reprex so that you can test the code. It should work...

If it doesn't work, here are some suggestions:

  1. Make sure all your libraries are up to date, especially tidyverse, mapview and sf. On my side, I run the code with the following versions: tidyverse 1.3.1, mapview 2.10.0 and sf 1.0.6
  2. Close all open documents in your R session and close R. Reopen a new R session and open only the file that contains the code to test.
  3. Load only the libraries needed to run the code.

Hope this helps. I'm crossing my fingers that these suggestions will unstuck you.

Reprex

  • Your data
library(tidyverse)
library(mapview)
library(sf)

lat <- c(300786, 301348,  301467, 302384, 304644, 304747, 304863, 304882)
lon <- c(4215918, 4215650, 4215784, 4216077, 4211303, 4211336, 4211395, 4211457)
geodeticDA <- c("NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83", "NAD83")
utmZone <- c("11N", "11N", "11N", "11N", "11N", "11N", "11N", "11N")
ID <- c("A", "A", "A", "A", "B", "B", "B", "B")

df <- data.frame(lat, lon, geodeticDA, utmZone, ID)


plot_local1 <- st_as_sf(df, 
                        coords = c("lat", "lon"),
                        crs = "+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ")

mapview(plot_local1, map.types = "Esri.NatGeoWorldMap") 

  • The code that should work...
polygon <- plot_local1 %>% 
  dplyr::group_by(ID) %>% 
  dplyr::summarize() %>% 
  sf::st_cast("POLYGON")

mapview(polygon, map.types = "Esri.NatGeoWorldMap") 

Created on 2022-03-05 by the reprex package (v2.0.1)

lovalery
  • 4,524
  • 3
  • 14
  • 28
  • 1
    It's finally working! Apparently I was working in an older version of R which wasn't compatible with the updated packages. I thought I had all the packages updated since R studio didn't list any updates but that wasn't the case. Thank you for your help :) – JulietheFoodie Mar 07 '22 at 21:43
  • @JulietheFoodie thank you very much for your feedback. Glad that I could help. Cheers. – lovalery Mar 07 '22 at 22:49