Как связать внешнюю область полигонов Вороного и пересечься с данными карты

Я пытаюсь визуализировать результаты процедуры кластеризации kmeans для следующих данных, используя voronoi polygons на карте США.

Создание следующих

Вопрос

Я хотел бы иметь возможность привязывать внешнюю область полигонов и пересекать полученную область с моей картой Соединенных Штатов, чтобы полигоны полностью представляли собой территорию США. Я не смог понять, как это сделать. Любая помощь приветствуется.

Ответы

Ответ 1

Моя конечная цель в задании этого вопроса заключалась в том, чтобы написать script, где я могу произвольно изменить число кластеров kmeans и быстро визуализировать результаты с помощью полигонов voronoi, которые покрывают мою желаемую область области.

Я еще не достиг этого, но я достиг достаточного прогресса, и решил, что опубликовать то, что у меня может привести к более быстрому решению.

# Create Input Data.Frame
input <- as.data.frame(cbind(x$long, x$lat))
colnames(input) <- c("long", "lat")

# Set Seed and Run Clustering Procedure
set.seed(123)
km <- kmeans(input, 35)

# Format Output for Plotting
centers <- as.data.frame(cbind(km$centers[,1], km$centers[,2]))
colnames(centers) <- c("long", "lat")
cent.id <- cbind(ID = 1:dim(centers)[1], centers)

# Create Spatial Points Data Frame for Calculating Voronoi Polygons
coords  <- centers[,1:2]
vor_pts <- SpatialPointsDataFrame(coords, centers, proj4string = CRS("+proj=longlat +datum=WGS84"))

Я также нашел ниже. функция во время поиска решения в Интернете.

# Function to Extract Voronoi Polygons 

SPdf_to_vpoly <- function(sp) {

  # tile.list extracts the polygon data from the deldir computation
  vor_desc <- tile.list(deldir([email protected][,1], [email protected][,2]))

  lapply(1:length(vor_desc), function(i) {

    # tile.list gets us the points for the polygons but we 
    # still have to close them, hence the need for the rbind

    tmp <- cbind(vor_desc[[i]]$x, vor_desc[[i]]$y)
    tmp <- rbind(tmp, tmp[1,])

    # Now we can make the polygons
    Polygons(list(Polygon(tmp)), ID = i)
  }) -> vor_polygons
  # Hopefully the caller passed in good metadata
  sp_dat <- [email protected]

  # This way the IDs should match up with the data & voronoi polys
  rownames(sp_dat) <- sapply(slot(SpatialPolygons(vor_polygons), 'polygons'), slot, 'ID')

  SpatialPolygonsDataFrame(SpatialPolygons(vor_polygons), data = sp_dat)
}

С помощью указанной функции определенные полигоны могут быть извлечены соответственно

vor     <- SPdf_to_vpoly(vor_pts)
vor_df  <- fortify(vor)

Чтобы получить полигоны voronoi, чтобы они хорошо соответствовали карте США, я загрузил cb_2014_us_state_20m с веб-сайта Census и побежал следующее:

# US Map Plot to Intersect with Voronoi Polygons - download from census link and place in working directory
us.shp <- readOGR(dsn = ".", layer = "cb_2014_us_state_20m")
state.abb <- state.abb[!state.abb %in% c("HI", "AK")]

Low48 <- us.shp[[email protected]$STUSPS %in% state.abb,]

# Define Area Polygons and Projections and Calculate Intersection
Low48.poly <- as(Low48, "SpatialPolygons")
vor.poly   <- as(vor, "SpatialPolygons")

proj4string(vor.poly) <- proj4string(Low48.poly)
intersect  <- gIntersection(vor.poly, Low48.poly, byid = T)


# Convert to Data Frames to Plot with ggplot
Low48_df <- fortify(Low48.poly)
int_df   <- fortify(intersect)

Отсюда я мог визуализировать свои результаты с помощью ggplot, как и раньше:

# Plot Results
StateMap <- ggplot() + geom_polygon(data = Low48_df, aes(x = long, y = lat, group = group), col = "white")

StateMap +
  geom_polygon(data = int_df, aes(x = long, y = lat, group = group, fill = id), alpha = .4) +
  geom_point(data = input, aes(x = long, y = lat), col = factor(km$cluster)) +
  geom_label(data = centers, aes(x = long, y = lat, label = row.names(centers)), alpha =.2) +
  scale_fill_hue(guide = 'none') +
  coord_map("albers", lat0 = 30, lat1 = 40)

Сводка обновлений