rtidyversepolygonr-sfmapview

Create Polygons by Category in sf file R


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!


Solution

  • 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

    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") 
    

    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)