Python数据可视化 | 11、基于地理特征的美国农药使用分析

最新推荐文章于 2023-11-17 12:00:00 发布

AI阅读和图谱

最新推荐文章于 2023-11-17 12:00:00 发布

阅读量663

点赞数 2

分类专栏：代码实战 python 文章标签： python 数据可视化

本文链接：https://blog.csdn.net/qq_34740277/article/details/119921453

版权

python 同时被 2 个专栏收录

18 篇文章 0 订阅

订阅专栏

代码实战

14 篇文章 2 订阅

订阅专栏

数据背景

本案例数据取自著名的数据科学竞赛网站 Kaggle 上的美国农药使用数据集（Pesticide Use in Agriculture）。这套数据集很有趣，它记录了2014和2015年度每个各州各县的 423 种农药的使用情况。此外，该数据汇总的农药使用情况有两种估测指标来描述，分别是低估测值（LOW_ESTIMATE）和高估测值（HIGH_ESTIMATE）。以下是对这两个不同估测方法的官方解释：

Two different methods were used to estimate a range of pesticide use for all states except California. Both low and high estimate methods incorporated proprietary surveyed rates for United States Department of Agriculture Crop Reporting Districts, but the estimates differed in how they treated situations when a district was surveyed and pesticide use was not reported. Low estimates assumed zero use in the district for that pesticide; however, high estimates treated the unreported use of pesticides as missing data and estimated the pesticide usage from neighboring locations within the same region.

简而言之，美国农业部在收集国内各州县地区的农药使用情况时，部分地区的调查结果会显示“未上报”的情况。所以，“低估测值”方法在遇到类似情况时，会假定该地区没有使用某农药，在数据集中标记 NaN；“高估测值”方法会将未上报的农药地区标记为数据缺失，然后将估算邻近地区的农药使用情况作为该地区的估测值。

在本案例中，我们将会充分的体会和挖掘数据集的潜在信息。目标不仅是要对数据集做清理和可视化，还要学会在依据一个完整思路完成对数据集的探索与理解。此外，还会接触到如何对地理性特征做可视化呈现。

本案例会涉及的 Python 扩展库：numpy，pandas，matplotlib，seaborn，以及处理地理信息用到的 Python 扩展库：geopandas，contextily，rasterio。

确认数据

# 导入扩展库
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns  # seaborn-0.9.0
%matplotlib inline

# 加载数据集
data2014 = pd.read_csv('./2014.csv')
data2015 = pd.read_csv('./2015.csv')
states = pd.read_csv('./dictionary.csv')

# 合并数据集
data2014_states = data2014.merge(states, how='left', on=['STATE_CODE', 'COUNTY_CODE'])
data2015_states = data2015.merge(states, how='left', on=['STATE_CODE', 'COUNTY_CODE'])

在上面的代码中，我们依据两列特征 STATE_CODE 和 COUNTY_CODE 将变量 states 中的特征 COUNTY 和 STATE 合并到了 data2014 和 data2015 两个数据集中。
这是因为 STATE_CODE 和 COUNTY_CODE 两列特征是三个数据集的共有列，并且可以将其转化为对应的 COUNTY 和 STATE 的两列信息。

到此，我们可以查看一下目前数据集的基本状态：

data2014_states.head()

在这里插入图片描述

data2015_states.head()

在这里插入图片描述
从上面的两个表可以看到两个年度数据集的特征列是一致的，从左至右分别农药成分（COMPOUND）、年份（YEAR）、州代码（STATE_CODE）、县代码（COUNTY_CODE）、低估测值（LOW_ESTIMATE）、高估测值（HIGH_ESTIMATE）、县名称（COUNTRY）和州名称（STATE）。
通过调用 pd.DataFrame.info 函数，我们可以发现低估测值 LOW_ESTIMATE 在两个年度数据集中都存在缺失值。
正如我们在数据背景的介绍中提到的，这是由于部分州县并没有上报农药使用情况所造成的，于是与之相对应的高估测值 HIGH_ESTIMATE 就必然是通过临近地区的农药使用情况测算的结果。
由于低估测值特征的缺省值有着特定的物理含义，加之对应的缺失样本占比较大，所以我们并不打算对其进行样本剔除等其他常规数据清理操作，而是直接将其替换为值 -1 表征其含义（不用值 0 是因为其本身已有“已上报且该农药使用未0”的含义）。

# 填充缺失值
data2014_states['LOW_ESTIMATE'] = data2014_states.LOW_ESTIMATE.fillna(-1)
data2015_states['LOW_ESTIMATE'] = data2015_states.LOW_ESTIMATE.fillna(-1)

# 查验
print(data2014_states[data2014_states.LOW_ESTIMATE.isnull() == True].size)
print(data2015_states[data2015_states.LOW_ESTIMATE.isnull() == True].size)

数据统计与可视化分析

这一节，我们将开始对数据集进行统计化和可视化分析。

首先，我们要将列特征信息进行概念分类。根据我们要考察的研究目标（农药）而言，农药的成分和估计值可以化作一个类别，使我们分析的指标中必然要包含的特征信息。同时，其中的高低估计值可以作为一对对比类别进行比较分析，类似地，所有列指标依据两个不同年份的对应比较，也可以形成比较分析。剩下的州县及其对应代码作为地理型特征是非常特殊的一类，可以特别为此绘制地理可视化。

对列特征类别进行概念分类后，作为唯一的数值型特征的高低估计值，我们先要考察它在不同年份、各州县的分布情况。不过值得留意的是一定要尽可能的将信息充分的“暴露”出来，因为只有将细微的差异放大才能有新的发现和理解。比如下图：

sns.catplot(x="STATE", y="HIGH_ESTIMATE", kind="bar", data=data2014_states)
fig = plt.gcf()
fig.set_size_inches(20,5)

在这里插入图片描述
上图中，我们绘制的是各州的农药使用高估测值的平均值和95%置信区间。可以发现横坐标多数州的农药平均估测值是不超过4000的，而且在为数不少的州中体现出过于宽泛的置信区域。所以，为了挖掘农药使用量在较低水平的细微差异，同时对比高低估测值的细微差异，可以如下绘制箱型图：

temp_l = data2014_states[data2014_states.LOW_ESTIMATE != -1].copy()
temp_l.loc[:,'Estimate_method'] = 'LOW_ESTIMATE'
temp_l.loc[:,'estimate'] = temp_l.LOW_ESTIMATE
temp_h = data2014_states
temp_h.loc[:,'Estimate_method'] = 'HIGH_ESTIMATE'
temp_h.loc[:,'estimate'] = temp_l.HIGH_ESTIMATE
temp = pd.concat([temp_l, temp_h], axis=0).sort_values('STATE')

plt.figure(figsize=(30,8))
sns.boxplot(x="STATE", y="estimate", hue='Estimate_method', data=temp)
plt.yscale('log')
plt.xlabel("State")
plt.ylabel("Usage pesticide-use")
plt.title("Pesticide-use in all states in year 2014")
plt.grid()

在这里插入图片描述
上图是2014年度的农业使用情况。可以发现，不论是哪一种估测结果，不同州的农药使用量都有着明显的不同。其中，超过 $10^6$ 农业用量的州县就出现在 CA、FL、ID 和 WA 四个州中（留意纵坐标是用 log 函数约化过的），并且其他所有州县的农药用量都远少于这四个州县，其中普遍农药用量相对偏少（最大量小于 $10^4$ ）的州县只有来自 CT、NH、RI 和 VT 四个州中。尽管如此，每个州中也都有不使用任何农药的州县(可以很容易验证，所有州中都存在农药估测值为0的县)。

在各州的不同区县中，尽管农药用量的差异很大，相差5到6个数量级，但是绝大部分区县的农药用量范围基本都在 $10^3$ 以下，甚至个别州中绝大部分区县的农药用量都不超过 $10^2$ 。

虽然根据前一节的考察，高低估测方法并不意味着估测值绝对大小的差异，但是由上图的分布表现可以得知高估计值方法在各州上的分布总是比低估计值方法要偏高，这可以从箱型图的中位线或上四分卫线的对比可以得知。

同理，我们也可以照猫画虎的给出2015年度的可视化图像。结论大同小异。

temp_l = data2015_states[data2015_states.LOW_ESTIMATE != -1].copy()
temp_l.loc[:,'Estimate_method'] = 'LOW_ESTIMATE'
temp_l.loc[:,'estimate'] = temp_l.LOW_ESTIMATE
temp_h = data2015_states
temp_h.loc[:,'Estimate_method'] = 'HIGH_ESTIMATE'
temp_h.loc[:,'estimate'] = temp_l.HIGH_ESTIMATE
temp = pd.concat([temp_l, temp_h], axis=0).sort_values('STATE')

plt.figure(figsize=(30,8))
sns.boxplot(x="STATE", y="estimate", hue='Estimate_method', data=temp)
plt.yscale('log')
plt.xlabel("State")
plt.ylabel("Usage pesticide-use")
plt.title("Pesticide-use in all states in year 2015")
plt.grid()

在这里插入图片描述

或许，上两图中纵坐标经过对数化的箱型图让你难以体会农药使用量在各州的差异变化，那么我们可以在上两图的分析基础上，进一步对应的给出两个年度各州的农药平均使用量：

states2014mean = data2014_states[data2014_states.LOW_ESTIMATE != -1][['STATE', 'LOW_ESTIMATE', 'HIGH_ESTIMATE']].groupby('STATE').mean()
states2015mean = data2015_states[data2015_states.LOW_ESTIMATE != -1][['STATE', 'LOW_ESTIMATE', 'HIGH_ESTIMATE']].groupby('STATE').mean()


states2014mean.plot(kind='bar',grid=True,figsize=(20,5))
plt.xticks(rotation=0)
plt.title("Mean value of pesticide-use of all states in year 2014")

states2015mean.plot(kind='bar',grid=True,figsize=(20,5))
plt.xticks(rotation=0)
plt.title("Mean value of pesticide-use of all states in year 2015")

在这里插入图片描述

若考虑两个年份之间在各州农用使用量的变化，可以更明确的将时间特征信息带来的差异充分体现出来：

(states2015mean - states2014mean).fillna(0).plot(kind='bar',grid=True,figsize=(20,5))
plt.xticks(rotation=0)
plt.title("The Differences (changes) in all States between year 2014 and year 2015")
plt.show()

在这里插入图片描述
上图描述的从2014年到2015年，各州的农药平均使用量在两种估测方法上的变化情况。其中由于没有 CA 州在2015年的数据，所以图中并未绘制其变化情况。值得留意的是，两种估测量在个别州的变化差异上出现了低估测值大于高估测值的现象。

到此，上述的可视化分析还远不能够让我们满足，因为各州中不同区县的农药使用信息被直接压缩掉了，并且我们还没有考虑究竟是哪些农药成分的在其中产生影响，以及两个年份的农药使用成分变化差异。

我们把地理信息放到后面去强调，接下来是探讨究竟是哪些农药成分是美国人民最爱用的，以及不同年份里农药的使用情况。

num = 10
startangle = 90

# LOW_ESTIMATE
compound2014sum = data2014_states[data2014_states.LOW_ESTIMATE != -1].groupby(['COMPOUND']).sum()
compound2015sum = data2015_states[data2015_states.LOW_ESTIMATE != -1].groupby(['COMPOUND']).sum()
labels2014 = compound2014sum.LOW_ESTIMATE.sort_values(ascending = False).head(num).index.append(pd.Index(['Others']))
fracs2014 = compound2014sum.LOW_ESTIMATE.sort_values(ascending = False).head(num)
fracs2014 = fracs2014.append(pd.Series(compound2014sum.LOW_ESTIMATE.sort_values(ascending = False).tail(-num).sum()))
labels2015 = compound2015sum.LOW_ESTIMATE.sort_values(ascending = False).head(num).index.append(pd.Index(['Others']))
fracs2015 = compound2015sum.LOW_ESTIMATE.sort_values(ascending = False).head(num)
fracs2015 = fracs2015.append(pd.Series(compound2015sum.LOW_ESTIMATE.sort_values(ascending = False).tail(-num).sum()))
explode= (0.05,) + (0,)*(num)
color = sns.color_palette("Paired")

plt.figure(1, figsize=(12,12))
plt.subplot(2,2,1)
plt.pie(fracs2014, explode=explode, labels=labels2014,
                autopct='%1.1f%%', shadow=True, startangle=startangle, colors = color)
plt.title('Compound in all states in year 2014', bbox={'facecolor':'0.8', 'pad':5})
plt.ylabel('LOW_ESTIMATE')

plt.subplot(2,2,2)
plt.pie(fracs2015, explode=explode, labels=labels2015,
                autopct='%1.1f%%', shadow=True, startangle=startangle, colors = color)
plt.title('Compound in all states in year 2015', bbox={'facecolor':'0.8', 'pad':5})

# HIGH_ESTIMATE
compound2014sum = data2014_states[data2014_states.LOW_ESTIMATE != -1].groupby(['COMPOUND']).sum()
compound2015sum = data2015_states[data2015_states.LOW_ESTIMATE != -1].groupby(['COMPOUND']).sum()
labels2014 = compound2014sum.HIGH_ESTIMATE.sort_values(ascending = False).head(num).index.append(pd.Index(['Others']))
fracs2014 = compound2014sum.HIGH_ESTIMATE.sort_values(ascending = False).head(num)
fracs2014 = fracs2014.append(pd.Series(compound2014sum.LOW_ESTIMATE.sort_values(ascending = False).tail(-num).sum()))
labels2015 = compound2015sum.HIGH_ESTIMATE.sort_values(ascending = False).head(num).index.append(pd.Index(['Others']))
fracs2015 = compound2015sum.HIGH_ESTIMATE.sort_values(ascending = False).head(num)
fracs2015 = fracs2015.append(pd.Series(compound2015sum.HIGH_ESTIMATE.sort_values(ascending = False).tail(-num).sum()))

plt.subplot(2,2,3)
plt.pie(fracs2014, explode=explode, labels=labels2014,
                autopct='%1.1f%%', shadow=True, startangle=startangle, colors = color)
plt.ylabel('HIGH_ESTIMATE')

plt.subplot(2,2,4)
plt.pie(fracs2015, explode=explode, labels=labels2015,
                autopct='%1.1f%%', shadow=True, startangle=startangle, colors = color)
plt.show()

在这里插入图片描述

上图描述的从2014年到2015年，在各州县上的平均使用变化最大的前25种农药成分在两种估测方法上的变化情况。
并不意外的发现，Sulfuric Acid 农药成分的热度在新的2015年度中明显下降，这也就间接地解释了在2015年度中，为什么 Glyphosate 农药成分既是总用量最大，也是在各州县中最普遍的原因（从2014年的州县平均使用量第二名到2015年度的第一名）。

接下来，我们开始着重考虑地理特征列。因为上述的可视化仅仅是对地理信息（州）进行了罗列与排序，然而地理位置还有经纬度等隐含坐标需要我们将其挖掘出来，并进一步观察各种农药成分及其各年度的农药使用估测量上的地理分布情况。

我们开始使用 geopandas 和 contextily 包来更方便的处理和美化地理空间数据。（安装详情请参照官网，若出现载入包 import geopandas 报错，可参考：ISSUE 解决问题）美国的地理经纬数据从这里下载：US Census 或 github。

geopandas / rasterio / mapclassify / descartes 安装办法：

conda install geopandas
conda install rasterio
conda install mapclassify
conda install descartes

import geopandas as gp

us_geod = gp.GeoDataFrame.from_file('./cb_2014_us_county_500k/cb_2014_us_county_500k.shp', encoding='UTF-8')
us_geod['COUNTY_CODE'] = us_geod.COUNTYFP.map(lambda x: int(x))
us_geod['STATE_CODE'] = us_geod.STATEFP.map(lambda x: int(x))

# Check
print([code for code in data2014_states.COUNTY_CODE.unique() if code not in us_geod.COUNTYFP.map(lambda x: int(x))])
print([code for code in data2014_states.STATE_CODE.unique() if code not in us_geod.STATEFP.map(lambda x: int(x))])
print([code for code in data2015_states.COUNTY_CODE.unique() if code not in us_geod.COUNTYFP.map(lambda x: int(x))])
print([code for code in data2015_states.STATE_CODE.unique() if code not in us_geod.STATEFP.map(lambda x: int(x))])

data2014_geod = gp.GeoDataFrame(data2014_states)   # 转换格式
data2014_merge = gp.GeoDataFrame.merge(data2014_geod, us_geod, on = ['COUNTY_CODE','STATE_CODE'], how = 'left') # 合并
data2015_geod = gp.GeoDataFrame(data2015_states)   # 转换格式
data2015_merge = gp.GeoDataFrame.merge(data2015_geod, us_geod, on = ['COUNTY_CODE','STATE_CODE'], how = 'left') # 合并

print(sum(np.isnan(data2014_merge['HIGH_ESTIMATE'])))
print(sum(np.isnan(data2015_merge['HIGH_ESTIMATE'])))
data2014_merge.head()

在这里插入图片描述
我们将下载的美国州县地理数据分别和两个年度的农药使用数据进行了合并，并且确保了有农药数据的美国州县都是有对应的经纬地理信息，同时考虑了未上报农药使用情况的州县和验证了数据合并后的正确性。

首先，我们先绘制出2014-2015年的农药使用总量高估测值在各州县的分布：

# import contextily as ctx
import rasterio as rio
rtr = rio.open('america_6.tif')

def add_basemap(ax, rtr):
    # NOTE the transpose of the image data
    img = np.array([ band for band in rtr.read() ]).transpose(1, 2, 0)
    bb = rtr.bounds
    xmin, xmax, ymin, ymax = (bb.left, bb.right, bb.bottom, bb.top)#ax.axis()

    ax.imshow(img, extent=(xmin, xmax, ymin, ymax))
    ax.axis((xmin, xmax, ymin, ymax))
    ax.set_axis_off()

reported2014_sum = gp.GeoDataFrame(pd.merge(data2014_merge.groupby(['STATE_CODE','COUNTY_CODE']).sum().reset_index(), us_geod[['STATE_CODE','COUNTY_CODE', 'geometry']], on = ['STATE_CODE','COUNTY_CODE'], how = 'left'))
reported2014_sum.crs = {'init': 'epsg:4326'}
reported2015_sum = gp.GeoDataFrame(pd.merge(data2015_merge.groupby(['STATE_CODE','COUNTY_CODE']).sum().reset_index(), us_geod[['STATE_CODE','COUNTY_CODE', 'geometry']], on = ['STATE_CODE','COUNTY_CODE'], how = 'left'))
reported2015_sum.crs = {'init': 'epsg:4326'}


fig, ax1 = plt.subplots(1, figsize = (20,20))
add_basemap(ax1, rtr)
reported2014_sum.to_crs(rtr.crs).plot('HIGH_ESTIMATE',scheme='Quantiles', k=300, cmap='OrRd', ax=ax1, alpha=0.5, edgecolor='k')
ax1.set_title('Usage pesticide-use (High estimate), US, 2014 ', fontsize = 15)#设置图形标题
plt.show()

fig, ax2 = plt.subplots(1, figsize = (20,20))
add_basemap(ax2, rtr)
reported2015_sum.to_crs(rtr.crs).plot('HIGH_ESTIMATE',scheme='Quantiles', k=300, cmap='OrRd', ax=ax2, alpha=0.5, edgecolor='k')
ax2.set_title('Usage pesticide-use (High estimate), US, 2015', fontsize = 15)#设置图形标题
plt.show()

在这里插入图片描述

上图中农药总量在各区县地理位置分布上的总体表现。我们可以明显看到不同地区的农药使用量有明显的区域化特点。而且，通过对比美国的地形图，我们可以发现在普遍大量使用农药的地区，基本对应于美国的河流沿岸和平原地带等较容易种植作物的地区。

最后，我们要考察两个年份里不同区县对那种农药成分使用量最大：（高估算量）

#import contextily as ctx
import rasterio as rio
rtr = rio.open('america_6.tif')

def add_basemap(ax, rtr):
    # NOTE the transpose of the image data
    img = np.array([ band for band in rtr.read() ]).transpose(1, 2, 0)
    bb = rtr.bounds
    xmin, xmax, ymin, ymax = (bb.left, bb.right, bb.bottom, bb.top)#ax.axis()

    ax.imshow(img, extent=(xmin, xmax, ymin, ymax))
    ax.axis((xmin, xmax, ymin, ymax))
    ax.set_axis_off()


table = data2014_merge.groupby(['STATE_CODE','COUNTY_CODE','COMPOUND'])['HIGH_ESTIMATE'].mean()
table = table.unstack('COMPOUND').reset_index().fillna(0)
table['compound'] = [table.columns[index] for index in table.values.argmax(axis=1)]
reported2014_com = gp.GeoDataFrame(pd.merge(table[['STATE_CODE','COUNTY_CODE','compound']], us_geod[['STATE_CODE','COUNTY_CODE', 'geometry']], on = ['STATE_CODE','COUNTY_CODE'], how = 'left'))
reported2014_com['compound_simple'] = reported2014_com['compound'].replace(reported2014_com['compound'].value_counts().tail(-16).index, 'Others')
reported2014_com.crs = {'init': 'epsg:4326'}

table = data2015_merge.groupby(['STATE_CODE','COUNTY_CODE','COMPOUND'])['HIGH_ESTIMATE'].mean()
table = table.unstack('COMPOUND').reset_index().fillna(0)
table['compound'] = [table.columns[index] for index in table.values.argmax(axis=1)]
reported2015_com = gp.GeoDataFrame(pd.merge(table[['STATE_CODE','COUNTY_CODE','compound']], us_geod[['STATE_CODE','COUNTY_CODE', 'geometry']], on = ['STATE_CODE','COUNTY_CODE'], how = 'left'))
reported2015_com['compound_simple'] = reported2015_com['compound'].replace(reported2015_com['compound'].value_counts().tail(-16).index, 'Others')
reported2015_com.crs = {'init': 'epsg:4326'}

fig, ax1 = plt.subplots(1, figsize = (20,20))
add_basemap(ax1, rtr)
reported2014_com.to_crs(rtr.crs).plot('compound_simple',categorical=True, ax=ax1,cmap='prism', alpha=0.4, edgecolor='k',legend = True)
ax1.set_title('Most usage pesticide-use of compound (High estimate), US, 2014 ', fontsize = 15)#设置图形标题
plt.show()

fig, ax2 = plt.subplots(1, figsize = (20,20))
add_basemap(ax2, rtr)
reported2015_com.to_crs(rtr.crs).plot('compound_simple',categorical=True, ax=ax2,cmap='prism', alpha=0.4, edgecolor='k',legend = True)
ax2.set_title('Most usage pesticide-use of compound (High estimate), US, 2015', fontsize = 15)#设置图形标题
plt.show()

在这里插入图片描述

从上面两图可以得知，农药成分使用量也呈现地理区域化分布特点，可能是由于邻近的地区气候环境适合种植一定类别的农作物等原因。此外，农药成分 Glyphosate 不仅是2014年和2015年使用量最大的农药配方，而且明显看到在美国中部广大地区和东部部分地区风靡已久。

小结

通过本案例的学习，可以熟练使用 Python 的各类扩展库对数据集有效的进行分析和操作，尤为对数据集的探索性分析产生很大的启发。
此外，深入接触了数据的高级统计分析和高级可视化操作，尤其是对地理信息的可视化提供了诸多宝贵的思路和素材。

AI阅读和图谱

关注

2
点赞
踩
6

收藏

觉得还不错? 一键收藏
2
评论
Python数据可视化 | 11、基于地理特征的美国农药使用分析

数据背景本案例数据取自著名的数据科学竞赛网站 Kaggle 上的美国农药使用数据集（Pesticide Use in Agriculture）。这套数据集很有趣，它记录了2014和2015年度每个各州各县的 423 种农药的使用情况。此外，该数据汇总的农药使用情况有两种估测指标来描述，分别是低估测值（LOW_ESTIMATE）和高估测值（HIGH_ESTIMATE）。以下是对这两个不同估测方法的官方解释：Two different methods were used to estimate a r
复制链接

扫一扫

专栏目录