如何使用R、matlab、arcgis进行趋势面分析

 

趋势面分析主要是用来研究变量在空间区域内的分布趋势。原理很简单，就是通过拟合模型z ~ f(x, y)（ f(x, y)一般指x和y的多项表达式 ），模拟地理要素的空间分布规律。然后通过得到的拟合结果，获得空间内任意点的预测值（可以理解为空间插值）

使用 ArcGis 分析

ArcGis本身提供两种类型的趋势分析方法：「趋势分析」和「趋势面分析」。

趋势分析

趋势分析工具位置如图

如果工具栏没有这个选项，则需要自己进行相关设置。具体步骤：菜单栏 > 自定义 > 扩展模块 > 勾选 「Geostatistical Analyst」 > 点击 「关闭」，然后鼠标右击工具栏空白处，勾选 「Geostatistical Analyst」。到此，就可以使用该工具进行相关分析。

进行趋势分析，如果是面数据，需要将其转换为点数据。工具位置：工具箱 > 系统工具箱 > 数据管理工具箱 > 要素 > 要素转点，具体操作步骤不再演示，见谅！

下面的趋势分析图形操作通过下面的视频进行演示（懒得写文字了😝）

最终结果如图

图中，两个曲线分别代表地理要素z与经度、纬度的拟合曲线（这里阶数设置为2），此外YOZ面、XOZ面、XOY面分别绘制了空间点的投影。

趋势面分析

通过搜索窗口快速定位趋势面分析工具，步骤如下

运行趋势面法工具，「输入点要素」选择我们需要分析的点文件，「z值字段」选择我们要分析的地理要素，「输出栅格」选择路径并命名文件名称。「多项式阶数」这里选为 2 阶，「回归类型」选择线性（若 z 为二值变量，可选LOGISTIC），点击确认。步骤及结果如图所示。

 

为了更加直观，自定义颜色分级显示，结果如图：

 

 

可以看出，该地理要素在空间分布上的规律还是比较明显的。趋势面分析与趋势分析原理其实是一样的，只不过趋势分析通过z ~ f(x)和z ~ f(y)得到拟合曲线，趋势面分析则通过z ~ f(x, y)得到拟合曲面，两者都可以反应地理要素的空间分布规律。

使用 matlab 分析

ArcGis虽然可以进行趋势面分析，但是它通过将得到的拟合曲面投影到平面上来反映空间分布规律。matlab则可以将拟合曲面在三维空间内进行展示，理解理解起来更加直观。

将各点的经纬度以及变量数据以列向量的形式导入matlab，下面直接给出计算过程：

% 绘制三维散点图
scatter3(lon, lan, fdi,'filled')
hold on;
% 估计模型，得到估计参数
lon2 = lon.^2;
lan2 = lan.^2;
lonlan = lon.*lan;
X = [ones(length(lon),1) lon lan lon2 lonlan lan2];
[b,bint,r,rint,stats] = regress(fdi,X,95);
% 生成格点经纬度，计算预测值
xfit = min(lon):0.1:max(lon);
yfit = min(lan):0.1:max(lan);
[XFIT,YFIT]= meshgrid (xfit,yfit);
ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT  + b(4) * XFIT.^2 + b(5) * XFIT .* YFIT + b(6) * YFIT.^2;
% 绘制拟合曲面
mesh(XFIT,YFIT,ZFIT);

最终图形如下：

使用 R 分析

这里主要使用R的plt3D程序包，下面给出代码及注释

setwd("E:/微信/3D图示")
library(readxl)
library(plot3D)

# 导入数据（包括经纬度以及地理要素数据）
data <- read_xls("经纬度.xls")
# 定义绘图函数
scatter3D_fancy <- function(x, y, z,..., colvar = z)
{
  # Compute the linear regression (z = ax + by + d)
  fit <- lm(z ~ x + y)
  # predict values on regular xy grid
  grid.lines = 40
  x.pred <- seq(min(x), max(x), length.out = grid.lines)
  y.pred <- seq(min(y), max(y), length.out = grid.lines)
  xy <- expand.grid( x = x.pred, y = y.pred)
  z.pred <- matrix(predict(fit, newdata = xy), 
                   nrow = grid.lines, ncol = grid.lines)
  # fitted points for droplines to surface
  fitpoints <- predict(fit)
  
  panelfirst <- function(pmat) {
    
    # 投影点 XOY
    XY <- trans3D(x, y, z = rep(min(z), length(z)), pmat = pmat)
    scatter2D(XY$x, XY$y, colvar = colvar, col = "brown", pch = 20, 
              cex = 1, add = TRUE, colkey = FALSE)
    
    # 投影点 XOZ
    XY <- trans3D(x, y = rep(min(y), length(y)), z, pmat = pmat)
    scatter2D(XY$x, XY$y, colvar = colvar, col = "purple", pch = 20, 
              cex = 1, add = TRUE, colkey = FALSE)
    
    # 投影点 YOZ
    XY <- trans3D(x = rep(max(x), length(x)), y, z, pmat = pmat)
    scatter2D(XY$x, XY$y, colvar = colvar, col = "orange", pch = 20, 
              cex = 1, add = TRUE, colkey = FALSE)
  }
  
  # 绘制主图
  scatter3D(x, y, z, ..., colvar = colvar, panel.first=panelfirst,
            colkey = list(length = 1, width = 0.7, cex.clab = 0.3, dist = -0.15),
            surf = list(x = x.pred, y = y.pred, z = z.pred,  
                        facets = NA, colvar = z.pred, fit = NULL), main = "趋势面分析") 
  
}
scatter3D_fancy(data1$lon, data1$lan, data1$外商直接投, phi = 20, bty = "g",  type = "h",
                ticktype = "detailed", theta = 200, d = 10, pch = ".", cex = 1)

最终结果如图：

 

更过精彩内容，点击关注微信公众号