1. 介紹
空間權(quán)重是許多空間分析的基礎(chǔ)膨处,可以用多種方式來(lái)表示。
PySAL能夠創(chuàng)建砂竖、處理和分析三類空間權(quán)重矩陣:
- Contiguity Based Weights
- Distance Based Weights
- Kernel Weights
2. PySAL空間權(quán)重類型
由于空間矩陣往往十分稀疏真椿,為了節(jié)省存儲(chǔ)空間,PySAL使用兩個(gè)字典來(lái)存儲(chǔ)空間矩陣——鄰居矩陣和權(quán)重矩陣乎澄。
2.1 鄰接權(quán)重
首先讓我們來(lái)創(chuàng)建一個(gè)5×5的網(wǎng)格(25個(gè)空間單位):
>>> import pysal
>>> w = pysal.lat2W(5, 5)
其中突硝,w對(duì)象有很多屬性。
>>> 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是空間單元的數(shù)目置济。
pct_nonzero表示矩陣的稀疏程度解恰。
weight用來(lái)存儲(chǔ)權(quán)重。
neighbor用來(lái)存儲(chǔ)鄰居浙于。
histogram則展示了鄰接情況护盈,上面的例子中,表示4個(gè)單元有2個(gè)鄰居羞酗,12個(gè)單元有3個(gè)鄰居腐宋,9個(gè)單元有4個(gè)鄰居。
鄰居的定義也是可以更改的。
>>> wq = pysal.lat2W(rook = False)
>>> wq.neighbors[0]
[5, 1, 6]
PySAL里有三種規(guī)則:Rook(車)脏款、Queen(后)围苫、Bishop(象)。
會(huì)玩國(guó)際象棋的話應(yīng)該比較容易理解撤师,Rook規(guī)則表示邊與邊相鄰才算鄰居剂府,Queen規(guī)則表示只要有一個(gè)角相鄰就算鄰居,Bishop規(guī)則表示有角相鄰且沒(méi)有邊相鄰才算鄰居剃盾。
Bishiop規(guī)則不是基本規(guī)則腺占,要通過(guò)PySAL的Difference方法計(jì)算得到。
lat2W
函數(shù)在需要使用規(guī)則形狀的矩陣來(lái)仿真時(shí)挺有用的痒谴。在實(shí)證研究中衰伯,通常會(huì)有一個(gè)shapefile,是一種非拓?fù)潢P(guān)系的向量數(shù)據(jù)結(jié)構(gòu)积蔚,然后就需要進(jìn)行和空間權(quán)重相關(guān)的空間分析意鲸。因?yàn)槲募胁缓負(fù)潢P(guān)系,因此在分析前尽爆,需要先建立空間權(quán)重矩陣怎顾。
PySAL會(huì)默認(rèn)根據(jù)鄰接圖建立權(quán)重,一般會(huì)用到.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規(guī)則漱贱,只要把上面代碼中的Rook換成Queen即可槐雾。
此外,還可以用含geometry列的dataframe建立權(quán)重矩陣幅狮。其中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)
再或者崇摄,也可以從可迭代的形狀對(duì)象直接建立擎值。
>>> 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
方法。但是它返回的對(duì)象是W
配猫,不是前面的Queen或者Rook幅恋。因?yàn)槿绻麤](méi)有原始形狀,不能確認(rèn)權(quán)重就是鄰接權(quán)重泵肄。
>>> W = pysal.weights.Rook.from_file(pysal.examples.get_path('columbus.gal'))
>>> type(W)
pysal.weights.weights.W
2.2 距離權(quán)重
除了鄰接,距離也可以用來(lái)確定權(quán)重淑翼。
注意:距離要在平面上計(jì)算腐巢,所以必須事先進(jìn)行投影。
2.2.1 k最近鄰權(quán)重
這里的鄰居用k最近鄰準(zhǔn)則定義玄括。
下面的例子中冯丙,用前面的5×5網(wǎng)格的質(zhì)心來(lái)作為測(cè)距的基點(diǎn)。首先,建兩個(gè)25維數(shù)組來(lái)存坐標(biāo)胃惜,放進(jìn)data里泞莉。
>>> import numpy as np
>>> x,y=np.indices((5,5))
>>> x.shape=(25,1)
>>> y.shape=(25,1)
>>> data=np.hstack([x,y])
接下來(lái)定義kNN權(quán)重:
>>> wknn3 = pysal.weights.KNN(data, k = 3)
>>> wknn3.neighbors[0]
[1, 5, 6]
>>> wknn3.s0
75.0
最近鄰計(jì)算用KD樹實(shí)現(xiàn)。為了避免多次重建KD樹船殉,提供了讓用戶改編權(quán)重的簡(jiǎn)便方法鲫趁。下面是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是最近鄰個(gè)數(shù)利虫,p是p范數(shù)挨厚。默認(rèn)是不新建kNN對(duì)象的。
其實(shí)我們也可以直接從shapefile建立kNN權(quán)重矩陣糠惫。
>>> 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 距離閾值權(quán)重
k-NN權(quán)重保證所有對(duì)象都有相同數(shù)量的鄰居疫剃。
另一種基于距離的權(quán)重則依靠距離閾值或距離段來(lái)定義,在研究對(duì)象周圍一定距離閾值內(nèi)的才是它的鄰居硼讽。
>>> 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]
可見巢价,不同對(duì)象可能鄰居不一樣多。
和前面一樣固阁,這個(gè)也能直接從dataframe創(chuàng)建:
>>> 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:(先確認(rèn)一下最小閾值)
>>> 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
距離閾值矩陣也可以讀取連續(xù)數(shù)值蹄溉。權(quán)重是距離的倒數(shù)。
>>> 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]
如果把衰減指數(shù)設(shè)為-2您炉,那結(jié)果就是重力權(quán)重:
>>> wid2 = pysal.weights.DistanceBand.from_array(points,14.2,alpha = -2.0, binary=False)
>>> wid2.weights[0]
[0.01, 0.0079999999999999984]
3. 核密度權(quán)重
將基于距離的權(quán)重和連續(xù)值權(quán)重結(jié)合起來(lái)柒爵,我們就得到了核密度權(quán)重。
>>> 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
參數(shù)相當(dāng)于距離的一個(gè)閾值赚爵,而核密度函數(shù)決定了衰減方式棉胀。 可選的核密度函數(shù)包括:‘triangular’, ’uniform’, ’quadratic’, ’epanechnikov’, ’quartic’, ’bisquare’, ’gaussian’。
bandwidth
是可以自定義的(可以是定值冀膝,也可以是可變的數(shù)組):
>>> 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]])
試試換一個(gè)核密度函數(shù):
>>> 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]])
具體內(nèi)容參見Kernel
唁奢。
類似地,也有Kernel.from_shapefile
和Kernel.from_dataframe
窝剖。