这个博客如果只是一个复刻 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$ ，被称之为“截断型”权重核函数。
空间距离
空间权重是根据距离计算的，而如何定义距离，也是个非常重要的问题。 我们最常用的就是欧氏距离（投影坐标系）或大圆距离（地理坐标系）。 除此之外，还有其他一些会用到的距离计算方法。
地理加权回归分析中，也经常用到以下几个距离：
- 曼哈顿距离：$D_{1,2}=|x_1-x_2|+|y_1-y_2|$
- 闵可夫斯基距离：$D_{1,2}=\left(|x_1-x_2|^p+|y_1-y_2|^p\right)^{\frac{1}{p}}$
- 路网距离：道路网络上两点的最短距离
对于栅格数据，或者栅格采样的点数据，也可以采用“四邻域距离”、“皇后距离”等等。
时空地理加权回归分析中，采用了一个“时空距离”的概念，即将时间和空间组合到一起。 原作者提供的思路是 $$ 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
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;
};
|
然后分别对 QgsdkWeight 和 QgsdkDistance 进行具体实现，如带宽权重和欧式/大圆距离（隐去 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 在使用的时候，就需要预先设置 mFocusPoints 和 mFocusPoints 属性。 如果指定 mGeographic 为 true 则按照大圆距离计算，反之按照欧氏距离计算。 这个 QgsdkBandwidthWeight 可以根据指定的不同的核函数进行计算。
