[PySAL] 空间权重

原文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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,496评论 6 501
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,407评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,632评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,180评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,198评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,165评论 1 299
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,052评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,910评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,324评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,542评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,711评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,424评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,017评论 3 326
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,668评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,823评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,722评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,611评论 2 353

推荐阅读更多精彩内容