此線程的擴展:Create choropleth map from coordinate points。 (爲了與儘可能多的人相關,我不想將這兩個線程結合起來。)R - Cloropleth:在多邊形內的數據點中,有多少百分比具有特定的列值?
我有一個由多個觀測值組成的數據幀,每個觀測值都有地理座標(緯度 - 經度)和布爾值(是 - 否)值。我想要生成一個世界的世界地圖,其中每個區域/多邊形都被其內相關布爾值等於true的點的百分比所着色。
這裏是一個最小可重現的例子,它現在只根據多邊形中的點數進行着色。 「喜歡」的數據列是我的布爾值。
# Load package
library(tidyverse)
library(ggmap)
library(maps)
library(maptools)
library(sf)
data <- data.frame(class = c("Private", "Private", "Private", "Private", "Private", "Private", "Not Private", "Not Private", "Private", "Private", "Not Private", "Private", "Not Private", "Private", "Private", "Not Private", "Not Private", "Private", "Private", "Not Private"),
lat = c(33.663944, 41.117936, 28.049601, 39.994684, 36.786042, 12.797659, 52.923318, 33.385555, 9.295242, 32.678207, 41.833585, -28.762956, 39.284713, 36.060964, 36.052239, 36.841764, 33.714237, 33.552863, 32.678207, -38.132401),
lon = c(-83.98686, -77.60468, -81.97271, -82.98577, -119.78246, 121.82814, -1.16057, -86.76009, 123.27758, -83.17387, -87.73201, 32.05737, -76.62048, -115.13517, -79.39961, -76.35592, -85.85172, -112.12468, -83.17387, 144.36946))
# Convert to simple feature object
point_sf <- st_as_sf(data, coords = c("lon", "lat"), crs = 4326)
# Get world map data
worldmap <- maps::map("world", fill = TRUE, plot = FALSE)
# Convert world to sp class
IDs <- sapply(strsplit(worldmap$names, ":"), "[", 1L)
world_sp <- map2SpatialPolygons(worldmap, IDs = IDs,
proj4string = CRS("+proj=longlat +datum=WGS84"))
# Convert world_sp to simple feature object
world_sf <- st_as_sf(world_sp)
# Add country ID
world_sf <- world_sf %>%
mutate(region = map_chr(1:length([email protected]), function(i){
[email protected][[i]]@ID
}))
# Use st_within
result <- st_within(point_sf, world_sf, sparse = FALSE)
# Calculate the total count of each polygon
# Store the result as a new column "Count" in world_sf
world_sf <- world_sf %>%
mutate(Count = apply(result, 2, sum))
# Convert world_sf to a data frame world_df
world_df <- world_sf
st_geometry(world_df) <- NULL
# Get world data frame
world_data <- map_data("world")
# Merge world_data and world_df
world_data2 <- world_data %>%
left_join(world_df, by = c("region"))
ggplot() +
geom_polygon(data = world_data2, aes(x = long, y = lat, group = group, fill = Count)) +
coord_fixed(1.3)
特別感謝https://stackoverflow.com/users/7669809/ycw尋求幫助。
謝謝你指出這些,我糾正了這個例子。 「雖然座標是經度/緯度,但假定它們是平面的」警告是沒有意義的,可以看着過去。 –