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