python 中提供了 KMeans库，可以方便我们对数据进行相应的聚类分析。
下面举个对于气温数据进行聚类分析的例子，数据来自ERA-5，可以自行从官网下载。
数据内容如下所示：
1、聚类分析
首先是导入库：
from sklearn.cluster import KMeans
然后对数据进行一下处理：
1、转换一下数据维度顺序，将高度level放到第一维
2、对nan值数据进行掩膜，
3、经纬度、时间、高度的大小，构造聚类矩阵
import xarray as xr import numpy as np path='D://t_data.nc' ds=xr.open_dataset(path).sel(lat=slice(-20,20)) da=ds.transpose('level','time','lat','lon') #转换一下数据维度顺序便于后续处理 time=da.time.data #读取时间 lat=da.lat.data #读取纬度 lon=da.lon.data #读取经度 omega=da.t.data #读取气温 w= omega.reshape(27,17*144*480) #构建矩阵，27表示level层数，其他几项代表经纬度、时间的size mask = np.isfinite(w[0,:]) #构建掩膜mask，判断是否有nan值 # # # # 依次循环将每一层得数据填入到零矩阵 # # 空矩阵k，便于后续聚类 k = np.zeros(((w[0,:][mask].shape[-1]),27)) # # 对于每一个高度层进行循环，填入每一层的非nan值数据 for i in range(27): k[:,i] = w[i,mask]
上面把矩阵构建好了，是个二维矩阵k，第2维表示高度层（level），第1维表示该层上所有的非nan值的数据，下面就可进行聚类了。K-means则是事先给定聚类数，这里的聚类数可以通过一些方法（SSE、轮廓系数）进行评估，以便选取适合的聚类数。
km=KMeans(n_clusters=6) # 6表示聚类的个数，聚类容器 km.fit(k) #进行聚类 label=km.labels_ #获取聚类的标签
等代码运行结束就可以发现聚类成功了，数据被分为0、1、2、3、4、5六类数据，这里可以将聚类后的数据单独保存，以便后续数据处理，直接导入即可，不需要再聚类一次浪费时间。（ps：这里提醒一下，对于同一套数据进行多次聚类，会导致label的标签序号每次可能都不同，但是个数是不变的，所以建议聚类一次后没问题就直接保存！！
2、保存数据
保存数据，可以直接使用下面两种方法都行，保存为.npy或者.nc，看个人喜好。
### 直接保存为.npy格式 np.save("/label.npy",label) ### 保存为 .nc格式 # label = time*lat*lon label_save = np.full((12*17*144), np.nan) label_save[mask] = label label_save = label_save.reshape(12,17,144) # ============================================================================================= # save labels label_nc = xr.Dataset( { "label":(("month","lat","lon"), label_save) }, coords={ "month":time, "lat":lat, "lon":lon, } ) label_nc.attrs["long_name"] = "K-mean labels" label_nc.to_netcdf("/Kmean6label.nc")
3、对于聚类的数据进行计算
聚类结束后，可以对于每一类数据进行计算，求平均值、标准差、绘制曲线等等。
# select type of omega label0 = np.array(np.where(label==0)) label1 = np.array(np.where(label==1)) label2 = np.array(np.where(label==2)) label3 = np.array(np.where(label==3)) label4 = np.array(np.where(label==4)) label5 = np.array(np.where(label==5)) # numlabel0=len(label0[0,:]) numlabel1=len(label1[0,:]) numlabel2=len(label2[0,:]) numlabel3=len(label3[0,:]) numlabel4=len(label4[0,:]) numlabel5=len(label5[0,:]) # 将二维数据转为一维 w0 = np.array([k[i,:] for i in label0]).squeeze() w1 = np.array([k[i,:] for i in label1]).squeeze() w2 = np.array([k[i,:] for i in label2]).squeeze() w3 = np.array([k[i,:] for i in label3]).squeeze() w4 = np.array([k[i,:] for i in label4]).squeeze() w5 = np.array([k[i,:] for i in label5]).squeeze() # calculate mean w0_mean = np.nanmean(w0, axis=0) w1_mean = np.nanmean(w1, axis=0) w2_mean = np.nanmean(w2, axis=0) w3_mean = np.nanmean(w3, axis=0) w4_mean = np.nanmean(w4, axis=0) w5_mean = np.nanmean(w5, axis=0) # calculate std w0_std = np.std(w0, axis=0) w1_std = np.std(w1, axis=0) w2_std = np.std(w2, axis=0) w3_std = np.std(w3, axis=0) w4_std = np.std(w4, axis=0) w5_std = np.std(w5, axis=0)
4、聚类个数评估
主要有以下两种方法对于聚类个数的大小进行评估：
1、SSE（簇内误方差）–手肘法
- SSE参数的核心思想就是计算误方差和，SSE的值越小，证明聚类效果越好，随着聚类数k的增大，样本划分会更加精细，每个簇的聚合程度会逐渐提高，那么误差平方和SSE自然会逐渐变小。并且，当k小于真实聚类数时，由于k的增大会大幅增加每个簇的聚合程度，故SSE的下降幅度会很大，而当k到达真实聚类数时，再增加k所得到的聚合程度回报会迅速变小，所以SSE的下降幅度会骤减，然后随着k值的继续增大而趋于平缓，也就是说SSE和k的关系图是一个手肘的形状，而这个肘部对应的k值就是数据的真实聚类数。
2、轮廓系数S
- 局部轮廓系数最大(极大值)，则该值为最优选择。对于下面的代码中，也就是说当s最大值所对应的k值，即为聚类的个数。
# method 1 SSE = [] for i in range(1,11): km = KMeans(n_clusters=i) km.fit(k) #获取K-means算法的SSE SSE.append(km.inertia_) plt.plot(SSE) plt.show() # method 2 from sklearn.metrics import silhouette_score S = [] # 存放轮廓系数 for i in range(2,10): kmeans = KMeans(n_clusters=i) # 构造聚类器 kmeans.fit(k) S.append(silhouette_score(k,kmeans.labels_,metric='euclidean')) # plt.plot(S) # plt.show()
一个努力学习python的海洋 水平有限，欢迎指正！！！ 欢迎评论、收藏、点赞、转发、关注。 关注我不后悔，记录学习进步的过程~~
