Create polygons by color

Question

I have x and y coordinates (latitude and longitude) and a third variable that indicates what color each region of the graph should be painted. I would like to generate polygons and save a shape file dividing the regions by color as shown below. Any simple way to do this in R?

Data can be downloaded at https://1drv.ms/u/s!AvopivTctLy2me53Fvso9JxJoq1t7A?e=wc3Doo

# figure
library(plotly)
plot_ly(data = aux, x = ~x, y = ~y, color = ~cluster2)

Answer 1

The solution was to do some conversions: data.frame > raster > SpatialPolygonsDataFrame > sf

dat <- aux
dat$cluster2 <- mgsub::mgsub(dat$cluster2, c('k1','k2','k3'), c('1','2','3'))
dat$cluster2 <- as.numeric(dat$cluster2)

# Create a raster layer
library(raster)
r <- rasterFromXYZ(dat)

# create SpatialPolygonsDataFrame by cluster
pol1 <- rasterToPolygons(r, fun=function(x){x==1}, dissolve = TRUE)
pol2 <- rasterToPolygons(r, fun=function(x){x==2}, dissolve = TRUE)
pol3 <- rasterToPolygons(r, fun=function(x){x==3}, dissolve = TRUE)

#join SpatialPolygonsDataFrame
ab <- rbind(pol1,pol2,pol3)

# sf object, split in Polygons and save
library(sf)
abc <- st_as_sf(ab)
abc <- st_cast(abc, "POLYGON")
st_write(abc, "filename", driver="ESRI Shapefile", layer_options = "OVERWRITE=true", update = TRUE) 

Answer 2

You need to infer the polygons in somehow, because you start with just a bunch of points.

This is how to get a shapefile from the convex hulls of the points per cluster:

library(tidyverse)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

data <- read_delim("~/Downloads/data.csv", delim = ";")
#> 
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#>   x = col_double(),
#>   y = col_double(),
#>   cluster2 = col_character()
#> )
data
#> # A tibble: 609 x 3
#>         x       y cluster2
#>     <dbl>   <dbl> <chr>   
#>  1 524385 7567995 k2      
#>  2 524395 7567995 k2      
#>  3 524405 7567995 k2      
#>  4 524415 7567995 k2      
#>  5 524425 7567995 k2      
#>  6 524435 7567995 k2      
#>  7 524445 7567995 k2      
#>  8 524455 7567995 k2      
#>  9 524465 7567995 k1      
#> 10 524475 7567995 k1      
#> # … with 599 more rows

points <- data %>% st_as_sf(coords = c("x", "y"))
points
#> Simple feature collection with 609 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 524325 ymin: 7567635 xmax: 524655 ymax: 7567995
#> CRS:           NA
#> # A tibble: 609 x 2
#>    cluster2         geometry
#>    <chr>             <POINT>
#>  1 k2       (524385 7567995)
#>  2 k2       (524395 7567995)
#>  3 k2       (524405 7567995)
#>  4 k2       (524415 7567995)
#>  5 k2       (524425 7567995)
#>  6 k2       (524435 7567995)
#>  7 k2       (524445 7567995)
#>  8 k2       (524455 7567995)
#>  9 k1       (524465 7567995)
#> 10 k1       (524475 7567995)
#> # … with 599 more rows
  
hulls <-
  data %>%
  nest(-cluster2) %>%
  transmute(
    cluster2,
    geometry = data %>% map(~ {
      .x %>%
        dplyr::select(x,y) %>%
        as.matrix() %>%
        st_multipoint(dim = "XY")
      })
  ) %>%
  st_as_sf() %>%
  st_convex_hull()
#> Warning: All elements of `...` must be named.
#> Did you want `data = c(x, y)`?
hulls
#> Simple feature collection with 3 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 524325 ymin: 7567635 xmax: 524655 ymax: 7567995
#> CRS:           NA
#> # A tibble: 3 x 2
#>   cluster2                                                              geometry
#> * <chr>                                                                <POLYGON>
#> 1 k2       ((524455 7567765, 524345 7567805, 524345 7567825, 524385 7567995, 52…
#> 2 k1       ((524455 7567685, 524335 7567765, 524335 7567785, 524435 7567975, 52…
#> 3 k3       ((524485 7567635, 524325 7567665, 524325 7567735, 524385 7567965, 52…

tibble() %>%
  ggplot(aes(color = cluster2)) +
    geom_sf(data = hulls, fill = "transparent") +
    geom_sf(data = points)

hulls %>% st_write("hulls.shp")
#> Writing layer `hulls' to data source `hulls.shp' using driver `ESRI Shapefile'
#> Writing 3 features with 1 fields and geometry type Polygon.

^{Created on 2021-10-04 by the reprex package (v2.0.1)}