本文介绍了R - Cloropleth：在多边形内的数据点中，有多少比例具有特定的列值？的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
此线程的扩展：。 （我不想将这两个线程结合起来，以便与尽可能多的人相关）。
我有一个由许多观察值组成的数据帧，每个都有地理坐标（纬度 - 经度）和布尔（是 - 否）值。我想要生成一个世界地图，其中每个区域/多边形都被其内相关布尔值等于true的点的百分比所遮蔽。
这是一个最小可重现的例子，它现在只根据多边形中的点数进行遮蔽。数据的like列是我的布尔值。
＃载入程序包
library（tidyverse）
库（ggmap）
库（映射）
库（maptools）
库（sf）
数据< - data.frame（class = c（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 ））
＃转换为简单要素对象
point_sf< - st_as_sf（data，coords = c（lon，lat），crs = 4326）
＃获取世界地图数据
worldmap< - maps :: map（world，fill = TRUE，plot = FALSE）
＃将世界转换为sp class
IDs< - sapply（strsplit（worldmap $ （世界地图，ID = ID，
proj4string = CRS（+ proj = longlat + datum = WGS84））$（$） b
$ b＃将world_sp转换为简单要素对象
world_sf< - st_as_sf（world_sp）
＃添加国家ID
world_sf< - world_sf％> ％
mutate（region = map_chr（1：length（world_sp @ polygons），function（i）{
world_sp @ polygons [[i]] @ ID
）））
＃使用st_within
结果< - st_within（point_sf，world_sf，sparse = FALSE）
＃计算每个多边形的总数
＃将结果存储为world_sf中的新列Count
world_sf< - world_sf％>％
mutate（Count = apply （result，2，sum））
＃将world_sf转换为数据框world_df
world_df< - world_sf
st_geometry（world_df）< - NULL
＃获取世界数据框
world_data< - map_data（world）
＃合并world_data和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）
特别感谢获取帮助。
解决方案
我们可以首先计算多边形中有多少个点，过滤数据集中标记为 Private code> class 列，然后再计算多边形中的多少个点。我们稍后可以通过使用私人计数除以所有计数并以100％重复计算百分比来计算百分比。
关于 sf 对象的一个很好的功能是它也是一个数据框。因此，从 dplyr 包中管理数据框的操作（例如 filter ）也适用于 sf 对象。因此，我们可以使用 point_private_sf< - point_sf％>％过滤器（%% Private中的类％）等命令轻松过滤这些点。
＃加载包
库（tidyverse）
库（映射）
库（maptools）
库（sf）
###数据准备
数据< - data.frame（class = c（Private，Private，Private，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.05 737，-76.62048，-115.13517，-79.39961，-76.35592，-85.85172，-112.12468，-83.17387,144.36946））
＃转换为简单要素对象
point_sf< - st_as_sf（数据，coords = c（lon，lat），crs = 4326）
＃获取世界地图数据
worldmap< - maps :: map（world，fill = TRUE，plot = FALSE）
＃将世界转换为sp类
IDs< - sapply（strsplit（worldmap $ names，：），[，1L）
world_sp< - map2SpatialPolygons（worldmap，IDs = ID，
proj4string = CRS（+ proj = longlat + datum = WGS84））
＃将world_sp转换为简单要素对象
world_sf< - st_as_sf（world_sp）
＃添加国家ID
world_sf< - world_sf％>％
mutate（region = map_chr（1：length（ world_sp @ polygons），function（i）{
world_sp @ polygons [[i]] @ ID
）））
###使用st_within进行分析
＃对所有点使用st_within
result_all< - st_within（point_sf，world_sf，sparse = FALSE）
＃在类列中通过Private过滤点数
point_private_sf< - point_sf％>％filter（class％in％Private）
＃使用st_within对于私人点
result_private< - st_within（point_private_sf，world_sf，sparse = FALSE）
###计算每个多边形的总数
＃将结果存储为ew列world_sf中的Count_all和Count_private
world_sf< - world_sf％>％
mutate（Count_all = apply（result_all，2，sum），
Count_private = apply（result_private，2 ，sum））％％>％
＃计算百分比
mutate（Percent = ifelse（Count_all == 0，Count_all，Count_private / Count_all * 100））
## ＃绘制数据
＃将world_sf转换为数据框world_df
world_df< - world_sf
st_geometry（world_df）< - NULL
＃获取世界数据框
world_data< - map_data（world）
＃合并world_data和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 =百分比））+
coord_fixed（1.3）
An extension of this thread: Create choropleth map from coordinate points. (I didn't want to combine the two threads for the sake of being relevant to as many people as possible.)
I have a data frame consisting of many observations, each with geocoordinates (latitude-longitude) and a boolean (yes-no) value. I would like to generate a choropleth map of the world where each region/polygon is shaded by the percentage of points within it where the associated boolean value is equal to true.
Here is a minimally reproducible example, which right now only shades according to the number of points in a polygon. The "like" column of data is my boolean.
# 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(world_sp@polygons), function(i){
world_sp@polygons[[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)
Special thanks to https://stackoverflow.com/users/7669809/ycw for the help thus far.
解决方案
We can first count how many points in a polygon, filter the datasets for records labeled as Private in the class column, and then count how many points in the polygon again. We can later calculate the percentage by using Private count number dividing by all count number and multipling by 100%.
One nice feature about the sf object is it is also a data frame. So operations to manage a data frame, such as filter from the dplyr package, also work for the sf object. So we can use command like point_private_sf <- point_sf %>% filter(class %in% "Private") to filter the points easily.
# Load package
library(tidyverse)
library(maps)
library(maptools)
library(sf)
### Data Preparation
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(world_sp@polygons), function(i){
world_sp@polygons[[i]]@ID
}))
### Use st_within for the analysis
# Use st_within for all points
result_all <- st_within(point_sf, world_sf, sparse = FALSE)
# Filter the points by "Private" in the class column
point_private_sf <- point_sf %>% filter(class %in% "Private")
# Use st_within for private points
result_private <- st_within(point_private_sf, world_sf, sparse = FALSE)
### Calculate the total count of each polygon
# Store the result as ew columns "Count_all" and "Count_private" in world_sf
world_sf <- world_sf %>%
mutate(Count_all = apply(result_all, 2, sum),
Count_private = apply(result_private, 2, sum)) %>%
# Calculate the percentage
mutate(Percent = ifelse(Count_all == 0, Count_all, Count_private/Count_all * 100))
### Plot the data
# 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 = Percent)) +
coord_fixed(1.3)
这篇关于R - Cloropleth：在多边形内的数据点中，有多少比例具有特定的列值？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！