基于栅格数据计算重心转移模型——以NDVI重心转移为例

重心转移模型可以反映对象在空间演变的过程，使用ArcGIS空间统计工具箱（Spatial Statistics Tools）中的平均中心（Mean Center）工具可以计算矢量数据的重心，栅格数据如何计算重心转移模型呢？本篇教程以NDVI重心转移为例给大家介绍一下如何使用R语言计算栅格数据的重心转移模型。

重心计算公式如下：

式中，为植被覆盖重心的经纬度坐标，、为像元i的几何重心，为像元i的像元值。

有了上面的公式，就可以在R语言中自定义函数进行计算，下面是一个1-12月NDVI重心转移模型的计算实例，内容包括：

• 重心转移模型计算
• 重心转移模型可视化
• 重心转移距离计算
• 重心转移距离可视化

栅格数据创建重心转移模型

实验数据获取

实验数据来源于【数据共享】1982–2020年中国5 km分辨率逐月NDVI数据集，其他类似的多波段栅格时间序列数据均可。

重心计算

数据的预处理，从上面获取得实验数据中提取12个月份的数据作为重心转移模型计算的输入。

library(terra)
library(tidyverse)

cdrndvi = rast("./cdrNDVI/cdr_MVC_CN.tif")
names(cdrndvi) = seq(as.Date("1982-01-01"), as.Date("2020-12-31"), "month")
cdrndvi
plot(cdrndvi)
cdrndvi12 = cdrndvi[[1:12]]
names(cdrndvi12) = seq(1:12)
cdrndvi12
plot(cdrndvi12)

提取了12个波段的数据：

12个月份数据的预览，图层名称已经改为月份：

• 重心转移模型R语言计算要点：
• 自定义函数，根据重心转移模型的公式写一个R语言函数
• terra包的as.data.frame可以将栅格转为数据框，并添加像元坐标信息，也就是获取、、这三个参数
• 数据的变形，宽数据转长数据
• 按月份分组计算重心

#计算重心的公式
gravity_xy = function(p, x, y){
  gx = sum(p*x)/sum(p)
  gy = sum(p*y)/sum(p)
  gxy = tibble(gx, gy)
  gxy
}

#应用重心计算公式计算年际重心变化
cdrmonth_gxy <- terra::as.data.frame(cdrndvi12, xy = T)%>% 
  pivot_longer(cols = starts_with("1"): starts_with("12"), 
               names_to = "Year", 
               values_to = "value") %>%
  mutate(Year = as.integer(Year)) %>% 
  group_by(Year) %>%   #按时间分组计算重心
  summarize(gxy = list(gravity_xy(value, x, y))) %>%   #分组计算
  unnest(gxy) 

write.csv(cdrmonth_gxy, file = "./Gravity/cdrmonth_gxy.csv") #输出重心计算结果

重心变化绘图

R语言绘制重心转移图

下面绘图代码详解可以参考往期推文：使用tidyverse绘制箭头相连的散点图

library(ggrepel)
library(ggsci)
#绘制月份重心变化图
p1 = ggplot(cdrmonth_gxy, aes(x = gx, y = gy, label = Year))+
  geom_segment(aes(
    xend=c(tail(gx, n=-1), NA), 
    yend=c(tail(gy, n=-1), NA)
  ),
  arrow=arrow(length=unit(0.3,"cm")),
  color = "#f6975e"
  )+
  geom_point(color = "#33a02c") +
  geom_text_repel(data=cdrmonth_gxy) +
  labs(x = "经度", y = "纬度")+
  theme_bw()+
  scale_color_npg()+
  scale_fill_npg()
p1
ggsave("./Rplots/逐月NDVI重心变化.jpg", p1, width = 16, height = 12, dpi = 600, units = "cm")
plotly::ggplotly(p1)

绘图结果见文末图片左半

QGIS绘制重心转移图

QGIS也可以很方便的将计算结果绘制成图，叠加一些在线影像可以很方便的看出重心转移的具体位置。下面是一些关键步骤：

在线底图使用的星图地球的数据源，数据空间位置准确，方便使用，详见：QGIS又多了一个在线底图新选择，星图地球数据云测评，超好用！真香！

重心移动距离

#计算逐月重心转移距离
#计算折线长度
# 将gx, gy转换为地理坐标，构建sf对象
library(sf)
cdrmonth_gxy_sf <- cdrmonth_gxy %>%
  st_as_sf(coords = c("gx", "gy"), crs = 4326)

# 计算两点之间的距离
calculate_distance <- function(data) {
  length_m <- st_distance(data)             #计算距离，单位为米
  data$length_m <- length_m
  return(data)
}

#计算逐月重心转移距离，单位米
cdrmonth_gxy_m = cdrmonth_gxy_sf %>% 
  group_by(Year) %>%
  summarize() %>%
  arrange(Year) %>%
  calculate_distance()   #使用自定义函数分组计算年份间的距离

cdrmonth_gxy_m = cdrmonth_gxy_m$length_m %>% 
  as_tibble() %>%                #距离计算结果为40*40矩阵，年际距离变化为第二条对角线
  rowid_to_column() %>%
  pivot_longer(-rowid, names_to = "colid", values_to = "value") %>%
  mutate(colid = as.numeric(str_extract(colid, "\\d+"))) %>%
  filter(colid == rowid + 1) %>%
  pull(value) %>% 
  as_tibble() %>% 
  mutate(Month = 2:12) %>% 
  mutate(Month = factor(Month, levels = 2:12))

#绘制距离变化图
library(units)
p2 = ggplot(cdrmonth_gxy_m, aes(x = Month, y = value))+
  geom_bar(stat = "identity")+
  theme_bw()+
  scale_color_npg()+
  scale_fill_npg()+
  labs(x ="月", y ="距离")
p2

library(patchwork)
p1+p2
ggsave("./Rplots/逐月NDVI重心变化和距离.jpg", p1+p2, width = 16, height =7, units = "cm", dpi = 600)

参考文献

1. 赵卓文, 张连蓬, 李行, 等. 基于MOD13Q1数据的宁夏生长季植被动态监测[J]. 地理科学进展, 2017, 36(6): 741–752.
2. https://r-graph-gallery.com/218-basic-barplots-with-ggplot2.html#width
3. 【数据共享】1982–2020年中国5 km分辨率逐月NDVI数据集
4. 使用tidyverse绘制箭头相连的散点图
5. QGIS又多了一个在线底图新选择，星图地球数据云测评，超好用！真香！
6. 推荐一个野外考察神器！+QGIS野外考察路径制图