QGIS 二次开发笔记(3)——空间距离和空间权重

这个博客如果只是一个复刻 QGIS 的教程的话,就没有什么价值,大家只要照着 QGIS 中的复制粘贴就可以了。所以这篇博客先来介绍一些 QGIS 中没有的。 我们使用 QGIS 做二次开发的目的,无非就是在软件中集成我们研发的一些算法,尤其是空间算法,不论是针对矢量数据还是栅格数据(我们主要研究矢量数据)。 对于矢量数据的空间算法而言,空间距离和空间权重非常重要,因为其反映了地理学第一定律:

任何事物都是与其他事物相关的,只不过相近的事物关联更紧密。

我们可以看到,不论是在莫兰指数(Moran’s I)中,还是空间自回归模型(SAR)中,或者地理加权回归模型(GWR)中,空间权重都是非常重要的。

空间权重

空间权重的计算可以是多种多样的,除了有一条总的原则:权重随距离的增加而减小,即权重是距离的单调减函数。该函数可称为“空间权重核函数”,即 $w(d)$。 该函数还可以引入一些参数,如 $w(d;b)$ ,参数 $b$ 可以事先指定,或者根据优化算法进行优化。

在莫兰指数、空间自回归模型等算法中,常常使用不含参数的空间权重函数,如“平方反距离函数”(忽略 $d=0$ 的情况) $$ w = \frac{1}{d^2} $$ 或者“指数反距离函数” $$ w = e^{-d} $$ 由于权重只有相对大小有意义,因此这些函数可以不用像概率密度函数一样要求 $\int w \ \mathrm{d}d = 1$ 。

在核密度分析、地理加权回归分析等算法中,常常使用含有一个参数 $b$ 的空间权重函数,该参数被称之为“带宽”(bandwidth)。 例如常用的高斯权重核函数 $$ w(d;b) = exp\left{ -\frac{1}{2} \left(\frac{d}{b}\right)^2 \right} $$ 该函数在任何距离上都有一定的权重,即使权重非常小,被称之为“非截断型”权重核函数。 还有一种权重核函数,如“双平方权重核函数” $$ w(d;b) = \left{ \begin{matrix} \left( 1 - \left( \frac{d}{b} \right)^2 \right)^2 & d \leq b \ 0 & d> b \end{matrix} \right. $$ 即超过带宽范围的位置上权重均为 $0$ ,被称之为“截断型”权重核函数。

空间距离

空间权重是根据距离计算的,而如何定义距离,也是个非常重要的问题。 我们最常用的就是欧氏距离(投影坐标系)或大圆距离(地理坐标系)。 除此之外,还有其他一些会用到的距离计算方法。

地理加权回归分析中,也经常用到以下几个距离:

  1. 曼哈顿距离:$D_{1,2}=|x_1-x_2|+|y_1-y_2|$
  2. 闵可夫斯基距离:$D_{1,2}=\left(|x_1-x_2|^p+|y_1-y_2|^p\right)^{\frac{1}{p}}$
  3. 路网距离:道路网络上两点的最短距离

对于栅格数据,或者栅格采样的点数据,也可以采用“四邻域距离”、“皇后距离”等等。

时空地理加权回归分析中,采用了一个“时空距离”的概念,即将时间和空间组合到一起。 原作者提供的思路是 $$ D_{1,2} = \mu((x_1-x_2)^2+(y_1-y_2)^2)+\lambda(t_1-t_2)^2 $$ 其中 $\mu+\lambda=1$ 且 $\mu,\lambda>0$ ,并选择合适的值以平衡时间和空间因素。

如果数据是线数据,那又该如何定义距离呢?目前有几种定义 Flow 距离的方式,但是总体上不是很令人满意。 而距离又需要满足非负性、同一性、对称性和三角不等式,因此往往需要根据实际研究内容的特点设计距离。

空间权重和距离的程序实现

根据以上介绍可以发现,空间权重离不开距离,各自又都多种多样,而且独立于算法。 在面向对象语言中,我们可以使用继承和多态特性,实现对空间权重和距离的封装:

  1. 将“空间权重”、“空间距离”分别定义为基类,再从其中派生出各种具体的权重和距离。
  2. 根据需要,将权重和距离进行组合,提供统一接口计算权重。

空间权重的声明如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
class QgsdkWeight
{
public:
enum WeightType
{
BandwidthWeight
};
public:
QgsdkWeight() {}
virtual ~QgsdkWeight() {}
virtual QgsdkWeight* clone() = 0;
public:
// 求权重向量
virtual vec weight(vec dist) = 0;
};

空间距离的声明如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
class QgsdkDistance
{
public:
enum DistanceType
{
CRSDistance,
MinkwoskiDistance,
DMatDistance
};
public:
explicit QgsdkDistance(int total) : mTotal(total) {};
QgsdkDistance(const QgsdkDistance& d) { mTotal = d.mTotal; };
virtual ~QgsdkDistance() {};
virtual QgsdkDistance* clone() = 0;
virtual DistanceType type() = 0;
int total() const;
void setTotal(int total);

public:
// 计算距离的函数
virtual vec distance(int focus) = 0;
// 返回结果的元素个数
virtual int length() const = 0;
// 求最大距离
double maxDistance();
// 求最小距离
double minDistance();

protected:
int mTotal = 0;
};

这里只是用一个 focus 整型变量来指定计算当前数据中第 $i$ 个点到其他所有点的距离。 在该设计下,需要将所有用于计算的数据保存在 QgsdkDistance 实例中,这是为了避免不同派生类所需参数不同的问题。 但也可以采用另一种方法,即接受一个 void * 类型的参数,这个参数指向计算所需要的所有数据,即可实现不同类型参数的传递。 或者也可以使用 Qt 中特有的 QVariant 类型,以在派生类中实现不同类型参数的传递。

然后可以构建一个组合类,将权重和距离组合起来:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
class QgsdkSpatialWeight
{
public:
QgsdkSpatialWeight();
QgsdkSpatialWeight(QgsdkWeight* weight, QgsdkDistance* distance);
QgsdkSpatialWeight(const QgsdkSpatialWeight& spatialWeight);
~QgsdkSpatialWeight();

QgsdkWeight *weight() const;
void setWeight(QgsdkWeight *weight);
void setWeight(QgsdkWeight& weight);

QgsdkDistance *distance() const;
void setDistance(QgsdkDistance *distance);
void setDistance(QgsdkDistance& distance);

public:
QgsdkSpatialWeight& operator=(const QgsdkSpatialWeight& spatialWeight);
QgsdkSpatialWeight& operator=(const QgsdkSpatialWeight&& spatialWeight);

public:
virtual vec weightVector(int i);
virtual bool isValid();

private:
QgsdkWeight* mWeight = nullptr;
QgsdkDistance* mDistance = nullptr;
};

然后分别对 QgsdkWeightQgsdkDistance 进行具体实现,如带宽权重和欧式/大圆距离(隐去 get/set 函数):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
class QgsdkBandwidthWeight : public QgsdkWeight
{
public:
typedef double (*KernelFunction)(double, double);
static KernelFunction Kernel[];
static double GaussianKernelFunction(double dist, double bw);
static double ExponentialKernelFunction(double dist, double bw);
static double BisquareKernelFunction(double dist, double bw);
static double TricubeKernelFunction(double dist, double bw);
static double BoxcarKernelFunction(double dist, double bw);

public:
QgsdkBandwidthWeight();
QgsdkBandwidthWeight(double size, bool adaptive, KernelFunctionType kernel);
QgsdkBandwidthWeight(const QgsdkBandwidthWeight& bandwidthWeight);
QgsdkBandwidthWeight(const QgsdkBandwidthWeight* bandwidthWeight);

virtual QgsdkWeight * clone() override;

public:
virtual vec weight(vec dist) override;

private:
double mBandwidth;
bool mAdaptive;
KernelFunctionType mKernel;
};

class QgsdkCRSDistance : public QgsdkDistance
{
public:
static vec SpatialDistance(const rowvec& out_loc, const mat& in_locs);
static vec EuclideanDistance(const rowvec& out_loc, const mat& in_locs);
static double SpGcdist(double lon1, double lon2, double lat1, double lat2);

public:
explicit QgsdkCRSDistance(int total, bool isGeographic);
QgsdkCRSDistance(const QgsdkCRSDistance& distance);

virtual QgsdkDistance * clone() override;

DistanceType type() override { return DistanceType::CRSDistance; }

public:
virtual vec distance(int focus) override;
int length() const override;

protected:
bool mGeographic = false;
mat* mFocusPoints = nullptr;
mat* mDataPoints = nullptr;
};

按照当前设计,这个 QgsdkCRSDistance 在使用的时候,就需要预先设置 mFocusPointsmFocusPoints 属性。 如果指定 mGeographictrue 则按照大圆距离计算,反之按照欧氏距离计算。 这个 QgsdkBandwidthWeight 可以根据指定的不同的核函数进行计算。

感谢您的阅读,本文由 HPDell 的个人博客 版权所有。如若转载,请注明出处:HPDell 的个人博客(http://hpdell.github.io/编程/qgisdev3-spatialweight/
QGIS 二次开发笔记(2)——显示图层
使用 GitLab Runner 为 R 包配置持续集成服务