这一篇纯粹是我的原创,各位大神复制黏贴的时候能不能标明下!
1、问题由来
我要做核密度估算相关的东西,同时希望将估算的结果绘制成图。
我发现有一个便捷的东东叫做“seaborn ”的包,它是一个基于matplotlib专门用于绘图数据统计图的,对于注重数据分析本身,而希望减少绘图操作的人来说是个福音。在seaborn包里绘制核密度图就是一行代码的事。
import seaborn as sns
sns.displot(dfdata, x="lon",y='lat',kind="kde", bw_adjust=.8, fill=True)#核密度图
displot()为绘图函数,其中参数kind中的kde表示核密度分布,其它参数大家可以参考官网。
绘制的结果如下,总体分析结果是可以接受的:
但是对于我来说,我要做核密度聚类的话,通过上述代码并不能得到核密度估算值得数字列表,因此我要解决的问题就是,既能得到核密度估算值列表,同时尽可能的绘制出它的分布图。
2、完整代码
我的数据格式是转换后的dataframe,里面有N个cloumn,这是一个空间数据集,其中我要用的是经度列lon和纬度列lat。重点解决的问题有两个:
- 进行核密度估算
- 绘制散点数据空间图
为了讲述方便,下面就直接上完整代码,然后我们再逐行进行分析:
import matplotlib.pyplot as plt
import numpy as np
from sklearn.neighbors import KernelDensity
from scipy.interpolate import griddata
x=df['lon'] #69行数据
y=df['lat'] #69行数据
#堆叠数据,将x和y堆叠成坐标对([x,y]),结果为shape为69*2的narray
point=np.vstack([x,y]).T#沿着第一个轴堆叠数组。
#计算核密度,窗宽0.2
kde = KernelDensity(kernel='epanechnikov', bandwidth=1).fit(point)#核密度估计
s=kde.score_samples(point)#数据点核密度
#空间坐标网格化,结果为shape为69*69的narray矩阵数据
xx,yy=np.meshgrid(np.linspace(x.min(),x.max(),200),np.linspace(y.min(),y.max(),200))
#通过griddata函数对z进行插值,返回的结果为69*69的narray
zz=griddata(points=point,values=s,xi=(xx,yy),method="linear") #对纵坐标也就是核密度进行插值
#应用contourf绘制等高线图
a=np.linspace(np.nanmin(zz),np.nanmax(zz),10)
myc=plt.contourf(xx,yy,zz,levels=a,alpha=0.75,cmap='Blues')#cmap=plt.cm.hot)
cbar = plt.colorbar(myc)
plt.show()
3、核密度估算
进行核密度估算用到sklearn机器学习包里的KernelDensity函数。KernelDensity是一个核密度估计函数。kernel density estimation是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Rosenblatt (1955)和Emanuel Parzen(1962)提出,又名Parzen窗(Parzen window)。Ruppert和Cline基于数据集密度函数聚类算法提出修订的核密度估计方法。函数如下:
KernelDensity ( * ,bandwidth= 1.0 , algorithm = 'auto' , kernel = 'gaussian' , metric = 'euclidean' , atol = 0 , rtol = 0 , width_first = True , leaf_size = 40 , metric_params = None )
这里有几个重点参数:
- kernel:{‘gaussian’, ‘tophat’, ‘epanechnikov’, ‘exponential’, ‘linear’, ‘cosine’}, default=’gaussian’。它一些供选择的核函数,也就是不同的核密度估算方法,默认的是高斯算法。
- metric:str, default=’euclidean’,要使用的距离度量。请注意,并非所有指标都适用于所有算法。有关可用算法的描述,BallTree请参阅文档 。KDTree请注意,密度输出的归一化仅对欧几里得距离度量是正确的。默认为“欧几里得”。
- bandwidth:窗宽,或者带宽,默认为1,是核密度估算精细度的一个参数。
KernelDensity 有以下几个方法:
翻译过来就是下面
其中的X就是我们要参与核密度估算的输入值,在本案例中就是上述代码中的point,它可以是一个一维列表list,也可时多维列表。我们这里应用了numpy的vstack()函数对lon和lat两个等长度列表进行叠加成一个narray,并通过取秩.T,得到的结果就是如下的数据结构,也就是坐标对:
[[119.34 31.24]
[122.46 34.61]
[123.47 36.16]
[119.76 32.76]
......
[119.97 31.85]
[119.03 33.41]
[121.81 33.49]
[118.83 33.49]]
通过fit(X)和score_samples(X)方法就能够得到我们要计算核密度估算值。我们原始的数据是69对坐标值,计算的结果s就是一个长度69的list列表。
4、matplotlib绘制散点图等高线
关于应用matplotlib绘制等高线,我网上查阅了很多资料,包括matplotlib官网,是走了不少弯路。
matplotlib绘制等高线一般用contourf()或者contour()方法,前者填充颜色,后者只绘制等高线条。这个并没有什么问题。但是我在网上查了一天,发现大多例子是这样的:
注意上图红笔圈出的部分,这些例子毫无例外,它们的z坐标都是通过数学公式计算出来的,并不是原始的离散数据,我上官网查了以下也是这个例子。那么对于离散型的数据该怎么样处理呢?
其实核心就两点:
- 1、空间坐标lon和lat网格化
- 2、对核密度估算值进行插值计算
4.1空间坐标网格化的弯路
在这一部分我是走了不少的弯路,应用numpy的meshgrid()函数即可解决,我也翻遍了很多例子,期初我的代码是这样的。
xx,yy=np.meshgrid(x,y)
随意绘制的结果就是这样的。
定眼一看,这还是用contour()没有填色的等高线线条图,这叫个乱啊,这哪里是等高线图嘛,分明就是一堆乱蜘蛛网。这该如何是好?
重点来了!!!
我发现xx,yy=np.meshgrid(x,y)中,原来的x,y列表并不是有序排列的,也不是均匀分布的。抱着试一试的态度,我想既然要空间网格化,我把网格化改成有序的是不是就不乱了?于是我用了下面一行代码。
xx,yy=np.meshgrid(np.linspace(x.min(),x.max(),200),np.linspace(y.min(),y.max(),200))
意思是说:将x,y分别取最小、最大值,并等间距取200个数值,这样产生出的网格空间就有序了啊。改正以后的图如下。
这才是我们希望的等高线嘛,所以,各位同学在进行空间网格化的时候一定,一定要记住meshgrid(x,y)里的x,y一定不能是原始无序的数据,而要把它写成等间距有序列表。网格化的本质是对空间画N多的小间隔,为数据插值服务的嘛。
4.2核密度插值计算
插值计算用到scipy包中的griddata()函数。
scipy.interpolate.griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)
points:具有形状 (n, D) 的浮点数的二维 ndarray,或具有形状 (n,) 的一维 ndarray 的长度 D 元组。(也就是我们通过vstack堆叠出的原始坐标对列表)
values:原始数据中的Z坐标,可以是标高h或者我们我们上述计算得到的x,y坐标点的核密度估算值列表,长度和lon,lat一样。
xi,:浮点数或复数的值ndarray,形状(n,D),这个是重点,是我们上述用网格化处理后的空间点数据,插值就往这些点上插。
method:{'linear', 'nearest', 'cubic'}, 可选,三种插值方法,默认为linear,至于这三者之间的区别其实还挺大的,一般选用默认,具体区别可以参考下面连接https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
4.3绘制核密度等高线(分布图)
数据处理好后,一切就水到渠成。下面三行代码解决:
myc=plt.contourf(xx,yy,zz,levels=a,alpha=0.75,cmap='Blues')#cmap=plt.cm.hot)
cbar = plt.colorbar(myc)
plt.show()
第一行绘图,第二行绘制色块图例,第三行显示图片。OK
注意:contourf()中的数据均使我们经过网格化和插值计算以后的数据,它们均为narray数据,维度(N,M)此时应该是一样的。xx表示经度,yy表示纬度,zz表示点核密度估算值。alpha颜色透明度,cmap为填充的颜色。
成果图如下:
在此说明下,我不是专业的IT人员啊,只是工作需要自己摸索着做,有不对的地方还往指点。另外,写这篇文章确实花了一定精力,搬运的话定要标明来处或给个链接,大家互相尊重以下。