问题描述
我拥有超过180,000个数据点的全局数据。对于每对坐标，我都有一个值，很多时候，我在同一坐标上有多个值（请参见数据框df的示例）。我想使用六角形网格将此数据绘制到世界地图上，并且一直在挣扎。 statbinhex选项ggplot不允许我将网格设置为1000平方千米，并且它仅计算hexbin中的点数，而不是hexbin中所有值的平均值，因此我切换到其他选项。现在，我正在尝试使用“ sp”包中的spsample函数对世界地图进行网格化，但是我一直遇到错误。
数据：
df<-structure(list(Z = c(3.23,3.518,1.961,0.9845,0.9845),latitude = c(-5.333,-19.01134381,-16.81667,-20.578928,11.068441,42.65,59.6111,55.8498,58.6388,57.0064,56.0202,57.2875,59.5252,63.8363,59.6032,60.6532,63.7764,59.3615,62.6603,58.9813,58.9988,58.984,63.5407,62.3942,59.36,59.7953,68.3704,57.3549,59.6068,63.6589,59.169,59.7762,56.6949,56.2811,61.6237,56.3035,56.7949,56.6454,65.5021,59.8754,59.0856,55.7247,56.7308,59.5479,56.7237,56.7821,58.5819,59.5112,67.8864,59.0272,58.9797,60.2414,59.0464,59.0805,59.7875,55.6308,42.64,42.534,42.60166667,41.2874,65.256706,42.68333333,42.61138889,47.12,63.3,49.13547,66.287,66.336,66.468,66.697,66.968,67.076,67.566,67.668,67.679,67.939,68.033,68.455,68.501,68.576,68.881,68.992,69.117,69.141,69.203,69.406,69.426,69.458,69.512
),longitude = c(141.6,146.2214001,145.6333,145.483398,76.509568,77.47,13.9202,14.2217,16.0795,14.4578,14.6835,17.9708,16.2606,20.127,13.8554,15.7272,20.8167,13.4909,17.399,15.1313,15.0579,14.7382,19.7277,17.7196,13.4549,17.5859,18.7693,15.762,13.8814,20.3222,15.1416,18.3233,13.1492,16.0232,17.4425,14.7285,16.5662,12.7691,22.0001,18.0014,14.6461,14.1954,13.0661,17.5769,12.8976,12.8581,14.8691,16.883,22.2536,14.9963,15.0096,14.48,15.0569,14.9042,17.6261,13.5288,2.09,2.465,1.093611111,24.6651,31.904297,1.205833333,1.063888889,6.63555,-150.5,2.63457,36.865,36.014,35.334,34.347,29.208,41.126,33.391,33.617,33.654,32.91,34.921,35.344,28.733,29.408,33.02,29.031,36.062,29.242,33.455,30.21,31.057,31.5,30.464),country = c("New Guinea","Australia","India","Kyrgyzstan","Sweden","France","Spain","Greece","Russia","USA","Russia")),row.names = c(NA,-100L),class = c("tbl_df","tbl","data.frame"))
世界地图：
library(rnaturalearth)
world <- ne_countries(scale = "medium",returnclass = "sf")
Hexbinninng世界地图（这是我的代码不起作用并且无法执行hexbinning的地方，我也不知道如何确保hexbins的大小为1000平方公里）：
size <- 100
hex_points <- spsample(world,type = "hexagonal",cellsize = size)
hex_grid <- HexPoints2Spatialpolygons(hex_points,dx = size)
Error in (function (classes,fdef,mtable) :
unable to find an inherited method for function ‘spsample’ for signature ‘"sf"’
我的计划是创建一个多边形，该多边形是一个六边形或世界的六边形网格，每个六边形为1000 km平方，然后将我的数据点与多边形形状文件相交，然后绘制该六边形内所有点的平均值在世界各地。
有人知道怎么做吗？
解决方法
sf package应该可以满足您的大部分需求。
如果要使用六边形大小相等的六边形网格，则必须使用等面积投影。说明差异：
world_6933 <- st_transform(world,6933)
world_grid_6933 <- st_make_grid(world_6933,n = c(100,100),what = 'polygons',square = FALSE,flat_topped = TRUE) %>%
st_as_sf() %>%
mutate(area = st_area(.))
p2 <- ggplot() +
geom_sf(data = world_grid_6933,aes(fill = units::drop_units(area))) +
geom_sf(data = world_6933,fill = NA,color = 'white')
cowplot::plot_grid(p1,p2,nrow = 2,labels = c('unequal hexagons','equal hexagons'))
在顶部图中，赤道附近的六边形明显更大。 （以m ^ 2为单位）
确定网格和投影后，请根据数据创建一个sf对象，并将其转换为相同的CRS投影。
# df <- from data pasted in your question
st_as_sf(df,coords = c('longitude','latitude'))
st_crs(df_sf) <- 4326
df_sf <- st_transform(df_sf,6933) # 6933 is the same crs as the hex grid
head(df_sf)
Simple feature collection with 6 features and 2 fields
geometry type: POINT
dimension: XY
bbox: xmin: 13662460 ymin: -2570656 xmax: 14108360 ymax: -679392.6
projected CRS: WGS 84 / NSIDC EASE-Grid 2.0 Global
# A tibble: 6 x 3
Z country geometry
<dbl> <chr> <POINT [m]>
1 3.23 New Guinea (13662457 -679392.6)
2 3.52 Australia (14108359 -2382208)
3 3.52 Australia (14051615 -2115478)
4 3.52 Australia (14037152 -2570656)
5 3.52 Australia (14037152 -2570656)
6 3.52 Australia (14037152 -2570656)
现在他们已经准备好加入和操纵了：
oined_sf <- st_join(df_sf,world_grid_6933,join = st_within) %>%
mutate(hex = floor(as.numeric(rownames(.)))) %>%
select(-area)
head(joined_sf)
Soined_sf <- st_join(df_sf,join = st_within) %>%
mutate(hex = floor(as.numeric(rownames(.)))) %>%
select(-area)
head(joined_sf)
Simple feature collection with 6 features and 3 fields
geometry type: POINT
dimension: XY
bbox: xmin: 13662460 ymin: -2570656 xmax: 14108360 ymax: -679392.6
projected CRS: WGS 84 / NSIDC EASE-Grid 2.0 Global
Z country hex geometry
1 3.230 New Guinea 720 POINT (13662457 -679392.6)
2 3.518 Australia 485 POINT (14108359 -2382208)
3 3.518 Australia 509 POINT (14051615 -2115478)
4 3.518 Australia 462 POINT (14037152 -2570656)
5 3.518 Australia 462 POINT (14037152 -2570656)
6 3.518 Australia 462 POINT (14037152 -2570656)
上面的（joined_df）是所有点，原始df中的Z和国家，以及它在world_grid_6933中所属于的六角线号。
这应该为您提供一个加入的专栏：
world_grid_6933$hex <- seq.int(nrow(world_grid_6933)
# Sum of Z by hexagon
summary <- joined_sf %>%
st_drop_geometry() %>%
group_by(hex) %>%
summarise(sum_z = sum(Z))
left_join(world_grid_6933,summary,by = 'hex') %>%
filter(!is.na(sum_z)) %>%
st_crop(box) %>% # box is a bounding box I made to zoom in on a relevant area
ggplot() +
geom_sf(aes(fill = sum_z)) +
geom_sf(data = st_crop(world_6933,box),fill = NA) +
geom_sf(data = st_crop(df_sf,color = 'white'[![enter image description here][2]][2])+
scale_fill_viridis_c(option = 'A')
