1

此线程的扩展：从坐标点创建等值线图。（为了与尽可能多的人相关，我不想将这两个线程结合起来。）

我有一个由许多观察组成的数据框，每个观察都有地理坐标（纬度-经度）和布尔值（是-否）。我想生成一个世界的等值线图，其中每个区域/多边形都被其中相关的布尔值等于 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(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)

特别感谢https://stackoverflow.com/users/7669809/ycw迄今为止的帮助。

4

1 回答 1

2

我们可以先计算多边形中有多少个点，过滤数据集以查找标记为列Private中的记录class，然后再次计算多边形中有多少个点。我们稍后可以通过使用Private计数除以所有计数并乘以 100% 来计算百分比。

该对象的一个​​不错的功能sf是它也是一个数据框。因此，管理数据框的操作（例如filter来自dplyr包的操作）也适用于sf对象。所以我们可以使用命令 likepoint_private_sf <- point_sf %>% filter(class %in% "Private")轻松过滤点。

# 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)
于 2017-08-27T11:50:17.523 回答