美文网首页我爱编程
[PySAL] 空间权重

[PySAL] 空间权重

作者: 红薯铺子 | 来源:发表于2018-04-07 22:57 被阅读0次

    原文here

    1. 介绍

    空间权重是许多空间分析的基础,可以用多种方式来表示。
    PySAL能够创建、处理和分析三类空间权重矩阵:

    • Contiguity Based Weights
    • Distance Based Weights
    • Kernel Weights

    2. PySAL空间权重类型

    由于空间矩阵往往十分稀疏,为了节省存储空间,PySAL使用两个字典来存储空间矩阵——邻居矩阵和权重矩阵。

    2.1 邻接权重

    首先让我们来创建一个5×5的网格(25个空间单位):

    >>> import pysal
    >>> w = pysal.lat2W(5, 5)
    

    其中,w对象有很多属性。

    >>> w.n
    25
    >>> w.pct_nonzero
    12.8
    >>> w.weights[0]
    [1.0, 1.0]
    >>> w.neighbors[0]
    [5, 1]
    >>> w.neighbors[5]
    [0, 10, 6]
    >>> w.histogram
    [(2, 4), (3, 12), (4, 9)]
    

    n是空间单元的数目。
    pct_nonzero表示矩阵的稀疏程度。
    weight用来存储权重。
    neighbor用来存储邻居。
    histogram则展示了邻接情况,上面的例子中,表示4个单元有2个邻居,12个单元有3个邻居,9个单元有4个邻居。

    邻居的定义也是可以更改的。

    >>> wq = pysal.lat2W(rook = False)
    >>> wq.neighbors[0]
    [5, 1, 6]
    

    PySAL里有三种规则:Rook(车)、Queen(后)、Bishop(象)。
    会玩国际象棋的话应该比较容易理解,Rook规则表示边与边相邻才算邻居,Queen规则表示只要有一个角相邻就算邻居,Bishop规则表示有角相邻且没有边相邻才算邻居。
    Bishiop规则不是基本规则,要通过PySAL的Difference方法计算得到。

    lat2W函数在需要使用规则形状的矩阵来仿真时挺有用的。在实证研究中,通常会有一个shapefile,是一种非拓扑关系的向量数据结构,然后就需要进行和空间权重相关的空间分析。因为文件中不含拓扑关系,因此在分析前,需要先建立空间权重矩阵。

    PySAL会默认根据邻接图建立权重,一般会用到.from_shapefile方法:

    >>> w = pysal.weights.Rook.from_shapefile(pysal.examples.get_path("columbus.shp")
    >>> w.n
    49
    >>> print "%.4f"%w.pct_nonzero
    0.0833
    >>> w.histogram
    [(2, 7), (3, 10), (4, 17), (5, 8), (6, 3), (7, 3), (8, 0), (9, 1)]
    

    如果用Queen规则,只要把上面代码中的Rook换成Queen即可。

    此外,还可以用含geometry列的dataframe建立权重矩阵。其中dataframe可以是geopandas的,也可用PySAL pandas IO extension的pdio。

    >>> import geopandas as gpd
    >>> test = gpd.read_file(pysal.examples.get_path('south.shp'))
    >>> W = pysal.weights.Queen.from_dataframe(test)
    >>> Wrook = pysal.weights.Rook.from_dataframe(test, idVariable='CNTY_FIPS')
    
    >>> pdiodf = pysal.pdio.read_files(pysal.examples.get_path('south.shp'))
    >>> W = pysal.weights.Rook.from_dataframe(pdiodf)
    

    再或者,也可以从可迭代的形状对象直接建立。

    >>> import geopandas as gpd
    >>> shapelist = gpd.read_file(pysal.examples.get_path('columbus.shp')).geometry.tolist()
    >>> W = pysal.weights.Queen.from_iterable(shapelist)
    

    最后还有.from_file方法。但是它返回的对象是W,不是前面的Queen或者Rook。因为如果没有原始形状,不能确认权重就是邻接权重。

    >>> W = pysal.weights.Rook.from_file(pysal.examples.get_path('columbus.gal'))
    >>> type(W)
    pysal.weights.weights.W
    

    2.2 距离权重

    除了邻接,距离也可以用来确定权重。

    注意:距离要在平面上计算,所以必须事先进行投影。

    2.2.1 k最近邻权重

    这里的邻居用k最近邻准则定义。

    下面的例子中,用前面的5×5网格的质心来作为测距的基点。首先,建两个25维数组来存坐标,放进data里。

    >>> import numpy as np
    >>> x,y=np.indices((5,5))
    >>> x.shape=(25,1)
    >>> y.shape=(25,1)
    >>> data=np.hstack([x,y])
    

    接下来定义kNN权重:

    >>> wknn3 = pysal.weights.KNN(data, k = 3)
    >>> wknn3.neighbors[0]
    [1, 5, 6]
    >>> wknn3.s0
    75.0
    

    最近邻计算用KD树实现。为了避免多次重建KD树,提供了让用户改编权重的简便方法。下面是reweight方法:

    >>> w4 = wknn3.reweight(k=4, inplace=False)
    >>> w4.neighbors[0]
    [1,5,6,2]
    >>> l1norm = wknn3.reweight(p=1, inplace=False)
    >>> l1norm.neighbors[0]
    [1,5,2]
    >>> set(w4.neighbors[0]) == set([1, 5, 6, 2])
    True
    >>> w4.s0
    100.0
    >>> w4.weights[0]
    [1.0, 1.0, 1.0, 1.0]
    

    其中,k是最近邻个数,p是p范数。默认是不新建kNN对象的。

    其实我们也可以直接从shapefile建立kNN权重矩阵。

    >>> wknn5 = pysal.weights.KNN.from_shapefile(pysal.examples.get_path('columbus.shp'), k=5)
    >>> wknn5.neighbors[0]
    [2, 1, 3, 7, 4]
    

    或者从dataframe建立:

    >>> import geopandas as gpd
    >>> df = gpd.read_file(ps.examples.get_path('baltim.shp'))
    >>> k5 = pysal.weights.KNN.from_dataframe(df, k=5)
    

    2.2.2 距离阈值权重

    k-NN权重保证所有对象都有相同数量的邻居。
    另一种基于距离的权重则依靠距离阈值或距离段来定义,在研究对象周围一定距离阈值内的才是它的邻居。

    >>> wthresh = pysal.weights.DistanceBand.from_array(data, 2)
    >>> set(wthresh.neighbors[0]) == set([1, 2, 5, 6, 10])
    True
    >>> set(wthresh.neighbors[1]) == set( [0, 2, 5, 6, 7, 11, 3])
    True
    >>> wthresh.weights[0]
    [1, 1, 1, 1, 1]
    >>> wthresh.weights[1]
    [1, 1, 1, 1, 1, 1, 1]
    

    可见,不同对象可能邻居不一样多。

    和前面一样,这个也能直接从dataframe创建:

    >>> import geopandas as gpd
    >>> df = gpd.read_file(pysal.examples.get_path('baltim.shp'))
    >>> pysal.weights.DistanceBand.from_dataframe(df, threshold=6, binary=True)
    

    或者shapefile:(先确认一下最小阈值)

    >>> thresh = pysal.min_threshold_dist_from_shapefile(pysal.examples.get_path("columbus.shp")
    >>> thresh
    0.61886415807685413
    >>> // 下面就可以获得 distance band weights 了:
    >>> wt = pysal.weights.DistanceBand.from_shapefile(pysal.examples.get_path("columbus.shp", threshold=thresh, binary=True)
    >>> wt.min_neighbors
    1
    >>> wt.histogram
    [(1, 4), (2, 8), (3, 6), (4, 2), (5, 5), (6, 8), (7, 6), (8, 2), (9, 6), (10, 1), (11, 1)]
    >>> set(wt.neighbors[0]) == set([1,2])
    True
    >>> set(wt.neighbors[1]) == set([3,0])
    True
    

    距离阈值矩阵也可以读取连续数值。权重是距离的倒数。

    >>> points = [(10, 10), (20, 10), (40, 10), (15, 20), (30, 20), (30, 30)]
    >>> wid = pysal.weights.DistanceBand.from_array(points,14.2,binary=False)
    >>> wid.weights[0]
    [0.10000000000000001, 0.089442719099991588]
    

    如果把衰减指数设为-2,那结果就是重力权重:

    >>> wid2 = pysal.weights.DistanceBand.from_array(points,14.2,alpha = -2.0, binary=False)
    >>> wid2.weights[0]
    [0.01, 0.0079999999999999984]
    

    3. 核密度权重

    将基于距离的权重和连续值权重结合起来,我们就得到了核密度权重。

    >>> points = [(10, 10), (20, 10), (40, 10), (15, 20), (30, 20), (30, 30)]
    >>> kw = pysal.Kernel(points)
    >>> kw.weights[0]
    [1.0, 0.500000049999995, 0.4409830615267465]
    >>> kw.neighbors[0]
    [0, 1, 3]
    >>> kw.bandwidth
    array([[ 20.000002],
           [ 20.000002],
           [ 20.000002],
           [ 20.000002],
           [ 20.000002],
           [ 20.000002]])
    

    bandwidth参数相当于距离的一个阈值,而核密度函数决定了衰减方式。 可选的核密度函数包括:‘triangular’, ’uniform’, ’quadratic’, ’epanechnikov’, ’quartic’, ’bisquare’, ’gaussian’。
    bandwidth是可以自定义的(可以是定值,也可以是可变的数组):

    >>> bw = [25.0,15.0,25.0,16.0,14.5,25.0]
    >>> kwa = pysal.Kernel(points,bandwidth = bw)
    >>> kwa.weights[0]
    [1.0, 0.6, 0.552786404500042, 0.10557280900008403]
    >>> kwa.neighbors[0]
    [0, 1, 3, 4]
    >>> kwa.bandwidth
    array([[ 25. ],
           [ 15. ],
           [ 25. ],
           [ 16. ],
           [ 14.5],
           [ 25. ]])
    

    可变的带宽可以是自生的:

    >>> kwea = pysal.Kernel(points,fixed = False)
    >>> kwea.weights[0]
    [1.0, 0.10557289844279438, 9.99999900663795e-08]
    >>> kwea.neighbors[0]
    [0, 1, 3]
    >>> kwea.bandwidth
    array([[ 11.18034101],
           [ 11.18034101],
           [ 20.000002  ],
           [ 11.18034101],
           [ 14.14213704],
           [ 18.02775818]])
    

    试试换一个核密度函数:

    >>> kweag = pysal.Kernel(points,fixed = False,function = 'gaussian')
    >>> kweag.weights[0]
    [0.3989422804014327, 0.2674190291577696, 0.2419707487162134]
    >>> kweag.bandwidth
    array([[ 11.18034101],
           [ 11.18034101],
           [ 20.000002  ],
           [ 11.18034101],
           [ 14.14213704],
           [ 18.02775818]])
    

    具体内容参见Kernel
    类似地,也有Kernel.from_shapefileKernel.from_dataframe

    相关文章

      网友评论

        本文标题:[PySAL] 空间权重

        本文链接:https://www.haomeiwen.com/subject/zlvyhftx.html