大数据分析——土地养分分析

一、选题的背景

我国土地辽阔，自然环境复杂多样，农业生产的种植制度也千差万别。土壤的供肥能力、作物的生长特性千差万别。作物种植、施肥量与土壤基础养分含量息息相关。

二、大数据分析设计方案

基于k-means算法对典型作物种植与施肥分区进行研究，并结合多种方法和手段，提取土地数据，随后通过查询农业资料，百度搜索等手段收集该省份典型作物的生长对土壤基础养分含量的要求。

三、数据分析步骤

1.数据源

Osgeo.cn

2.数据清洗

import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 
data = pd.read_csv(r'shanxi.csv',encoding='utf-8',index_col=['地市名','县名'])
describe=data.describe()   
describe1=data.isnull().sum()   
describe2=data.dropna()  
describe3=data.drop_duplicates()   
print(data)

3.大数据分析过程及采用的算法

利用scipy和sklearn模块，k-means算法，处理数据

from numpy.random import uniform
from scipy.spatial.distance import cdist
data = pd.read_csv(r'shanxi.csv',encoding='utf-8',index_col=['地市名','县名'])
#数据去空
data=data.dropna()  
#霍普金斯统计量计算，input:DataFrame类型的二维数据，output:float类型的霍普金斯统计量
#默认从数据集中抽样的比例为0.3
def hopkins_statistic(data:pd.DataFrame,sampling_ratio:float = 0.3) -> float: 
  #抽样比例超过0.1到0.5区间任意一端则用端点值代替
  sampling_ratio = min(max(sampling_ratio,0.1),0.5)
  n_samples = int(data.shape[0] * sampling_ratio)
  sample_data = data.sample(n_samples) #原始数据抽样后剩余的数据
  data = data.drop(index = sample_data.index) #,inplace = True)
  #原始数据中抽取的样本与最近邻的距离之和
  data_dist = cdist(data,sample_data).min(axis = 0).sum()
  #人工生成的样本点，从平均分布中抽样(artificial generate samples)
  ags_data = pd.DataFrame({col:uniform(data[col].min(),data[col].max(),n_samples)\
  for col in data})
  #人工样本与最近邻的距离之和
  ags_dist = cdist(data,ags_data).min(axis = 0).sum()
  #计算霍普金斯统计量H
  H_value = ags_dist / (data_dist + ags_dist)
  return H_value
#代码实现
print(hopkins_statistic(data))  #结果大于0.5即存在簇结构

# 阈值法主观判断聚类的个数
from scipy.cluster.hierarchy import linkage,fclusterfrom sklearn.preprocessing import MinMaxScaler
matplotlib.rc("font",family='Microsoft YaHei',weight="bold")

data = pd.read_csv(r'shanxi.csv',encoding='UTF-8',index_col=['地市名','县名'])
mm = MinMaxScaler()
mm_data = mm.fit_transform(data)

plt.figure(figsize=(20,12),dpi=80)
z0=linkage(mm_data, method='average', metric='euclidean')
for n in range(5,13):
    f0=fcluster(z0, t=n, criterion='maxclust')
    f_0=f0.tolist()
    f_0.sort()
    count=[]
    for i in range(n):
        count.append(0)
    for j in f_0:
        for m in range(n):
            if j==m+1:
                count[m]=count[m]+1
    print(count)
    # x:count的长度（最后取最大即可） y:不同分类的每一类的县的个数
    x=[]
    for q in range(n):
        x.append(q+1)
    annotation = plt.annotate(('x='+str(x[n-1]), 'count='+str(count[n-1])), xy=(x[n-1]+0.1, count[n-1]-2),textcoords='data', horizontalalignment="left",
                                arrowprops=dict(arrowstyle="simple",connectionstyle="arc3,rad=-0.1", alpha=0.9),#绘制箭头
                                bbox=dict(boxstyle="round", facecolor="w",edgecolor="0.5", alpha=0.1))
      
    plt.plot(x, count,label="分为{0}类的情况".format(n))

plt.grid(alpha=0.4, linestyle=':')
plt.legend(loc="best")
plt.xticks(x,rotation=45)

plt.xlabel("分类")
plt.ylabel("每一分类，不同类的个数（从大到小）")
plt.title("山西省不同分类数的情况对比图")

plt.show()
#分析可得大概可分为7-11类

#进一步分析分类数:
from sklearn.cluster import KMeans
plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 
data = pd.read_csv(r'shanxi.csv',encoding='UTF-8',index_col=['地市名','县名'])
mean_distortions = []
K = range(7,12)
for k in K:
    k_means = KMeans(n_clusters=k, random_state=0,n_init=20,max_iter=500)
    k_means.fit(data)
    mean_distortions.append(sum(np.min(cdist(data, k_means.cluster_centers_, metric='euclidean'), axis=1))/ data.shape[0])
plt.plot(K, mean_distortions, 'bx-')
plt.xlabel('k')
plt.ylabel(u'平均畸变程度')
plt.title(u'用肘部法确定最佳的K值')
plt.show()
#由图可知k为9最好

# 聚类质量（轮廓系数）
from sklearn import metrics
score=[]
for k in K:
    k_means = KMeans(n_clusters=k, random_state=0,n_init=20,max_iter=500) 
    k_means.fit(data)
    labels = k_means.labels_
    score.append(metrics.silhouette_score(data, labels, metric='euclidean'))
plt.plot(K, score, 'bx-')
plt.xlabel('k')
plt.ylabel(u'轮廓系数')
plt.title(u'轮廓系数变化图')
plt.show()
#由图可知k为9最好

由上述图像可知，分为九类最好。

plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 

data = pd.read_csv(r'shanxi.csv',encoding='utf-8',index_col=['地市名','县名'])
mm = MinMaxScaler()
mm_data = mm.fit_transform(data)
#参数的设置由上可知,类个数为9
k_means = KMeans(n_clusters=9, random_state=0,n_init=20,max_iter=500)
answer=k_means.fit_predict(data)
for i in range(len(answer)):
    answer[i]=answer[i]+1
print(answer)  #分类结果

data1=pd.read_csv(r'shanxi.csv',encoding='utf-8')
answer_0=answer.tolist()
answer_0=pd.DataFrame(answer_0,columns=["分类"])
answer_0=data1.join(answer_0)
answer_0=answer_0.sort_values(by="分类",ascending=True)
answer_0['地市名-县名'] = answer_0['地市名'] + "-" + answer_0['县名']
answer_0.to_csv("分析后数据.csv",index=0)  #得到数据分析后的数据

4.数据可视化

data = pd.read_csv(r'分析后数据.csv',encoding='UTF-8',index_col = '地市名')
x=data.loc["临汾市"]
y1=x["有机质"]
y2=x["有效磷"]
y3=x["全氮"]
y4=x["速效钾"]
y5=x["PH"]
x11=x["县名"].tolist()

plt.figure(figsize=(20,8),dpi=80)
bar_width = 0.15
x1 = list(range(len(x)))
x2 =[i+bar_width for i in x1]
x3 =[i+bar_width*2 for i in x1]
x4 =[i+bar_width*3 for i in x1]
x5 =[i+bar_width*4 for i in x1]

plt.bar(x1,y1,width=bar_width,label='有机质')
plt.bar(x2,y2,width=bar_width,label='有效磷')
plt.bar(x3,y3,width=bar_width,label='全氮')
plt.bar(x4,y4,width=bar_width,label='速效钾')
plt.bar(x5,y5,width=bar_width,label='PH')

plt.legend() 
plt.xticks(x3,x11)

plt.xlabel("临汾市")
plt.ylabel("有机质，全氮，有效磷，速效钾，PH")
plt.title("山西省临汾市与土壤基础养分条形图")
plt.grid(alpha=0.4, linestyle=':')
plt.show()

由此图可以看出，在临汾市的所有市区中，侯马市和曲沃县土壤的速效钾和有机质含量都是相对较高的。但是总体来说，临汾市各个地区所有的土壤基础养分含量相差不大，分布比较均匀，且土壤速效钾含量均在150mg/kg以上。据此，我们可以根据农作物的生长习性去确定作物应该种植在哪一个县区中，结合农作物的生长习性可以得出临汾市更适合种植玉米、马铃薯和小麦。

pd_data=pd.read_csv(r'分析后数据.csv',encoding='UTF-8')
pd_data1=pd_data.loc[:,["有机质","全氮","有效磷","速效钾","PH"]]

plt.figure(dpi=80)
scatter_matrix(pd_data1,figsize=(15,12))
plt.show()

此图展示了有机质，全氮，有效磷，速效钾，PH这五种土壤基础养分互相之间的关系变化，采用散点图+矩阵图的形式，让我们可以很形象的展现了五种土壤基础养分之间的关系。

df=pd.read_csv(r'分析后数据.csv',encoding='UTF-8')
df=df.loc[:,["有机质","全氮","有效磷","速效钾","PH"]]

plt.figure(figsize=(5,7), dpi= 300)
font0 = FontProperties(family='SimSun')
sns.heatmap(df.corr(), xticklabels=df.corr().columns, yticklabels=df.corr().columns, cmap='RdYlGn', center=0, annot=True)

plt.title('土壤养分相关性图', fontsize=20,fontproperties=font0)
plt.xticks(fontsize=5,fontproperties=font0)
plt.yticks(fontsize=5,fontproperties=font0)
plt.show()

根据此图可以看出，有机质和全氮、有效磷，速效钾，全氮和有效磷、速效钾，有效磷和速效钾的相关系数绝对值都大于0.3且小于0.8，认为有弱的相关性。其余基础养分之间没有相关性。

data = pd.read_csv(r'分析后数据.csv',encoding='UTF-8')
y1=data["有机质"]
y2=data["全氮"]
y3=data["有效磷"]
y4=data["速效钾"]
y5=data["PH"]
x=data["地市名-县名"]

plt.figure(figsize=(25, 20), dpi=80)
plt.scatter(x,y1,label="有机质")
plt.scatter(x,y2,label="全氮")
plt.scatter(x,y3,label="有效磷")
plt.scatter(x,y4,label="速效钾")
plt.scatter(x,y5,label="PH")

plt.legend(loc="best")
plt.xticks(x,rotation=90)

plt.title('山西省地市名-县名与基础养分的散点图')
plt.xlabel('地市名-县名')
plt.ylabel('有机质，全氮，有效磷，速效钾，PH')

plt.show()

散点图中可以看出，对于相同分类里的各个县和区，它们的各土壤基础养分值相差不大，可以看出其中的相关性。而对于不同的分类，它们的值相差较大，散点比较杂乱，规律性较差，没有较大的关联性。通过分析不同地市的土壤基础养分情况，可以确定在这些地区中适合种植哪些作物。

5.附完整程序源代码

import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 
data = pd.read_csv(r'shanxi.csv',encoding='utf-8',index_col=['地市名','县名'])
describe=data.describe()   
describe1=data.isnull().sum()   
describe2=data.dropna()  
describe3=data.drop_duplicates()   
print(data)

from numpy.random import uniform
from scipy.spatial.distance import cdist
data = pd.read_csv(r'shanxi.csv',encoding='utf-8',index_col=['地市名','县名'])
#数据去空
data=data.dropna()  
#霍普金斯统计量计算，input:DataFrame类型的二维数据，output:float类型的霍普金斯统计量
#默认从数据集中抽样的比例为0.3
def hopkins_statistic(data:pd.DataFrame,sampling_ratio:float = 0.3) -> float: 
  #抽样比例超过0.1到0.5区间任意一端则用端点值代替
  sampling_ratio = min(max(sampling_ratio,0.1),0.5)
  n_samples = int(data.shape[0] * sampling_ratio)
  sample_data = data.sample(n_samples) #原始数据抽样后剩余的数据
  data = data.drop(index = sample_data.index) #,inplace = True)
  #原始数据中抽取的样本与最近邻的距离之和
  data_dist = cdist(data,sample_data).min(axis = 0).sum()
  #人工生成的样本点，从平均分布中抽样(artificial generate samples)
  ags_data = pd.DataFrame({col:uniform(data[col].min(),data[col].max(),n_samples)\
  for col in data})
  #人工样本与最近邻的距离之和
  ags_dist = cdist(data,ags_data).min(axis = 0).sum()
  #计算霍普金斯统计量H
  H_value = ags_dist / (data_dist + ags_dist)
  return H_value
#代码实现
print(hopkins_statistic(data))  #结果大于0.5即存在簇结构
　　

# 阈值法主观判断聚类的个数
from scipy.cluster.hierarchy import linkage,fcluster
from sklearn.preprocessing import MinMaxScaler
matplotlib.rc("font",family='Microsoft YaHei',weight="bold")

data = pd.read_csv(r'shanxi.csv',encoding='UTF-8',index_col=['地市名','县名'])
mm = MinMaxScaler()
mm_data = mm.fit_transform(data)

plt.figure(figsize=(20,12),dpi=80)
z0=linkage(mm_data, method='average', metric='euclidean')
for n in range(5,13):
    f0=fcluster(z0, t=n, criterion='maxclust')
    f_0=f0.tolist()
    f_0.sort()
    count=[]
    for i in range(n):
        count.append(0)
    for j in f_0:
        for m in range(n):
            if j==m+1:
                count[m]=count[m]+1
    print(count)
    # x:count的长度（最后取最大即可） y:不同分类的每一类的县的个数
    x=[]
    for q in range(n):
        x.append(q+1)
    annotation = plt.annotate(('x='+str(x[n-1]), 'count='+str(count[n-1])), xy=(x[n-1]+0.1, count[n-1]-2),textcoords='data', horizontalalignment="left",
                                arrowprops=dict(arrowstyle="simple",connectionstyle="arc3,rad=-0.1", alpha=0.9),#绘制箭头
                                bbox=dict(boxstyle="round", facecolor="w",edgecolor="0.5", alpha=0.1))
      
    plt.plot(x, count,label="分为{0}类的情况".format(n))

plt.grid(alpha=0.4, linestyle=':')
plt.legend(loc="best")
plt.xticks(x,rotation=45)

plt.xlabel("分类")
plt.ylabel("每一分类，不同类的个数（从大到小）")
plt.title("山西省不同分类数的情况对比图")

plt.show()
#分析可得大概可分为7-11类

#进一步分析分类数:
from sklearn.cluster import KMeans
plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 
data = pd.read_csv(r'shanxi.csv',encoding='UTF-8',index_col=['地市名','县名'])
mean_distortions = []
K = range(7,12)
for k in K:
    k_means = KMeans(n_clusters=k, random_state=0,n_init=20,max_iter=500)
    k_means.fit(data)
    mean_distortions.append(sum(np.min(cdist(data, k_means.cluster_centers_, metric='euclidean'), axis=1))/ data.shape[0])
plt.plot(K, mean_distortions, 'bx-')
plt.xlabel('k')
plt.ylabel(u'平均畸变程度')
plt.title(u'用肘部法确定最佳的K值')
plt.show()
#由图可知k为9最好

# 聚类质量（轮廓系数）
from sklearn import metrics
score=[]
for k in K:
    k_means = KMeans(n_clusters=k, random_state=0,n_init=20,max_iter=500) 
    k_means.fit(data)
    labels = k_means.labels_
    score.append(metrics.silhouette_score(data, labels, metric='euclidean'))
plt.plot(K, score, 'bx-')
plt.xlabel('k')
plt.ylabel(u'轮廓系数')
plt.title(u'轮廓系数变化图')
plt.show()
#由图可知k为9最好

plt.rcParams['font.sans-serif']=['SimHei'] 
plt.rcParams['axes.unicode_minus']=False 

data = pd.read_csv(r'shanxi.csv',encoding='utf-8',index_col=['地市名','县名'])
mm = MinMaxScaler()
mm_data = mm.fit_transform(data)
#参数的设置由上可知,类个数为9
k_means = KMeans(n_clusters=9, random_state=0,n_init=20,max_iter=500)
answer=k_means.fit_predict(data)
for i in range(len(answer)):
    answer[i]=answer[i]+1
print(answer)  #分类结果

data1=pd.read_csv(r'shanxi.csv',encoding='utf-8')
answer_0=answer.tolist()
answer_0=pd.DataFrame(answer_0,columns=["分类"])
answer_0=data1.join(answer_0)
answer_0=answer_0.sort_values(by="分类",ascending=True)
answer_0['地市名-县名'] = answer_0['地市名'] + "-" + answer_0['县名']
answer_0.to_csv("分析后数据.csv",index=0)  #得到数据分析后的数据

data = pd.read_csv(r'分析后数据.csv',encoding='UTF-8',index_col = '地市名')
x=data.loc["临汾市"]
y1=x["有机质"]
y2=x["有效磷"]
y3=x["全氮"]
y4=x["速效钾"]
y5=x["PH"]
x11=x["县名"].tolist()

plt.figure(figsize=(20,8),dpi=80)
bar_width = 0.15
x1 = list(range(len(x)))
x2 =[i+bar_width for i in x1]
x3 =[i+bar_width*2 for i in x1]
x4 =[i+bar_width*3 for i in x1]
x5 =[i+bar_width*4 for i in x1]

plt.bar(x1,y1,width=bar_width,label='有机质')
plt.bar(x2,y2,width=bar_width,label='有效磷')
plt.bar(x3,y3,width=bar_width,label='全氮')
plt.bar(x4,y4,width=bar_width,label='速效钾')
plt.bar(x5,y5,width=bar_width,label='PH')

plt.legend() 
plt.xticks(x3,x11)

plt.xlabel("临汾市")
plt.ylabel("有机质，全氮，有效磷，速效钾，PH")
plt.title("山西省临汾市与土壤基础养分条形图")
plt.grid(alpha=0.4, linestyle=':')
plt.show()

pd_data=pd.read_csv(r'分析后数据.csv',encoding='UTF-8')
pd_data1=pd_data.loc[:,["有机质","全氮","有效磷","速效钾","PH"]]

plt.figure(dpi=80)
scatter_matrix(pd_data1,figsize=(15,12))
plt.show()

df=pd.read_csv(r'分析后数据.csv',encoding='UTF-8')
df=df.loc[:,["有机质","全氮","有效磷","速效钾","PH"]]

plt.figure(figsize=(5,7), dpi= 300)
font0 = FontProperties(family='SimSun')
sns.heatmap(df.corr(), xticklabels=df.corr().columns, yticklabels=df.corr().columns, cmap='RdYlGn', center=0, annot=True)

plt.title('土壤养分相关性图', fontsize=20,fontproperties=font0)
plt.xticks(fontsize=5,fontproperties=font0)
plt.yticks(fontsize=5,fontproperties=font0)
plt.show()

data = pd.read_csv(r'分析后数据.csv',encoding='UTF-8')
y1=data["有机质"]
y2=data["全氮"]
y3=data["有效磷"]
y4=data["速效钾"]
y5=data["PH"]
x=data["地市名-县名"]

plt.figure(figsize=(25, 20), dpi=80)
plt.scatter(x,y1,label="有机质")
plt.scatter(x,y2,label="全氮")
plt.scatter(x,y3,label="有效磷")
plt.scatter(x,y4,label="速效钾")
plt.scatter(x,y5,label="PH")

plt.legend(loc="best")
plt.xticks(x,rotation=90)

plt.title('山西省地市名-县名与基础养分的散点图')
plt.xlabel('地市名-县名')
plt.ylabel('有机质，全氮，有效磷，速效钾，PH')

plt.show()