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

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

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

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

在莫兰指数、空间自回归模型等算法中，常常使用不含参数的空间权重函数，如“平方反距离函数”（忽略 d=0 的情况） w=1d2 或者“指数反距离函数” w=e−d 由于权重只有相对大小有意义，因此这些函数可以不用像概率密度函数一样要求 ∫w dd=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. 曼哈顿距离：D1,2=|x1−x2|+|y1−y2|
2. 闵可夫斯基距离：D1,2=(|x1−x2|p+|y1−y2|p)1p
3. 路网距离：道路网络上两点的最短距离

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

时空地理加权回归分析中，采用了一个“时空距离”的概念，即将时间和空间组合到一起。 原作者提供的思路是 D1,2=μ((x1−x2)2+(y1−y2)2)+λ(t1−t2)2 其中 μ+λ=1 且 μ,λ>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;
};

然后分别对 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 可以根据指定的不同的核函数进行计算。