基于Earth Mover's Distance计算.tif栅格空间利用相似度的问题

我正在用动态布朗桥运动模型分析声学接收阵列的动物追踪数据，数据集包含动物ID、检测时间戳、经纬度，以及studyperiod和sgroup两个分组变量：

Elasmo            Datetime      Lat       Lon studyperiod sgroup
1 X10141 2019-02-13 15:29:00 25.72441 -79.30922      IS2019  naive
2 X10141 2019-02-13 15:44:00 25.72441 -79.30922      IS2019  naive
3 X10141 2019-02-13 15:48:00 25.72441 -79.30922      IS2019  naive
4 X10141 2019-02-13 17:17:00 25.72441 -79.30922      IS2019  naive
5 X10141 2019-02-13 17:20:00 25.72441 -79.30922      IS2019  naive
6 X10141 2019-02-13 18:00:00 25.72441 -79.30922      IS2019  naive

我先用GitHub包dBBMMhomeRange计算个体利用分布（UD），再经缩放、加权求和得到组水平UD并保存为.ascii栅格。之后在ArcMap中统一栅格范围并导出为tif，导入R后的栅格属性如下：

r1 <- raster("~/N_IS2017_sc.tif")
r4 <- raster("~/N_IS2020_sc.tif")

> r1;r4
class      : RasterLayer 
dimensions : 1450, 1264, 1832800  (nrow, ncol, ncell)
resolution : 50, 50  (x, y)
extent     : -31658.67, 31541.33, -36085.92, 36414.08  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : N_IS2017_sc.tif 
names      : N_IS2017_sc 
values     : 0, 0.0003713508  (min, max)

class      : RasterLayer 
dimensions : 1450, 1264, 1832800  (nrow, ncol, ncell)
resolution : 50, 50  (x, y)
extent     : -31658.67, 31541.33, -36085.92, 36414.08  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : N_IS2020_sc.tif 
names      : N_IS2020_sc 
values     : 0, 0.0004588088  (min, max)

现在想用Kranstauber等人（2016）提出的EMD计算组间空间利用相似度，使用move包的emdDists()函数：

## 创建栅格栈
allrasters <- stack(r1,r4)

## 计算EMD
emdDists <- emd(allrasters/cellStats(allrasters, sum), threshold = 700)

但运行后R总是崩溃，有两个疑问：

1. 文献里用的是UDStack对象输入emdDists()，不确定上述tif栅格能不能用来计算EMD？
2. 如果方法可行，有哪些降低R计算资源消耗的办法？

问题1：tif栅格能否用于EMD计算 #

可以用，但需要满足几个核心前提：

• 你已经完成归一化：代码里allrasters/cellStats(allrasters, sum)的操作是对的，确保每个UD的总和为1，符合EMD对概率分布的要求
• 栅格空间属性完全匹配：你的r1和r4的行列数、分辨率、范围完全一致，这是EMD计算的关键——它要求两个分布在完全相同的网格上计算距离，这点你已经满足
• 虽然文献用UDStack，但move包的emd()本质可接受RasterStack输入，UDStack只是带元数据的特殊RasterStack，核心数据结构一致

注意：你的栅格crs为NA，建议先赋予正确坐标系——EMD计算的距离基于空间坐标，无坐标系可能导致距离逻辑出错，甚至引发崩溃。

问题2：降低计算资源消耗的方法 #

你的栅格有近200万单元格，EMD本质是求解最优运输问题，对内存和算力要求极高，试试这些方法：

• 降低栅格分辨率：用aggregate()合并单元格，比如把50m分辨率改成100m，大幅减少单元格数量：

  r1_agg <- aggregate(r1, fact=2, fun=sum)
  r4_agg <- aggregate(r4, fact=2, fun=sum)
  # 重新归一化
  allrasters_agg <- stack(r1_agg/cellStats(r1_agg, sum), r4_agg/cellStats(r4_agg, sum))
  

• 过滤低概率单元格：把UD中概率极低（比如小于1e-6）的单元格设为0，减少非零计算单元：

  allrasters_filtered <- clamp(allrasters, lower=1e-6, useValues=FALSE)
  allrasters_filtered <- allrasters_filtered/cellStats(allrasters_filtered, sum)
  

• 优化内存使用：
  • 用memory()查看当前内存限制，Windows下用memory.limit()提升可用内存，Linux/macOS可设置环境变量调整
  • 用rm()清理无用变量后运行gc()释放内存，避免同时加载不必要的对象
• 分批次计算：如果有多个栅格对比，不要一次性全栈输入，两两单独计算EMD
• 更换求解器：move包的emd()依赖transport包，可直接调用transport::emd()，尝试设置solver="networkflow"，可能比默认求解器更高效

内容的提问来源于stack exchange，提问作者smok