这里我们用FAO(Food and Agriculture Organization)组织提供的数据集，练习一下如何利用python进行探索性数据分析。

探索性数据分析(Exploratory Data Analysis，简称EDA)是指对已有的数据(特别是调查或观察得来的原始数据)在尽量少的先验假定下进行探索，通过作图、制表、方程拟合、计算特征量等手段探索数据的结构和规律的一种数据分析方法，特别是当我们面对大数据时代到来的时候，各种杂乱的“脏数据”，往往不知所措，不知道从哪里开始了解目前拿到手上的数据时候，探索性数据分析就非常有效。

1、导包

我们先导入需要用到的包

%matplotlib inline
%config InlineBackend.figure_format='retina'
import matplotlib as mpl
from matplotlib import pyplot as plt
import seaborn as sns
import pandas as pd
import os,sys
import warnings 
sns.set_context("poster",font_scale=1.3)

接下来，加载数据集

#导入数据
data = pd.read_csv('../数据集/农粮数据/aquastat.csv.gzip',compression='gzip')#compression直接使用磁盘上的压缩文件
data.head()

看一下数据量，

data.shape
#输出结果
(143280, 7)

看一下数据的信息，

data.info()
输出结果
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 143280 entries, 0 to 143279
Data columns (total 7 columns):
country          143280 non-null object
region           143280 non-null object
variable         143280 non-null object
variable_full    143280 non-null object
time_period      143280 non-null object
year_measured    96411 non-null float64
value            96411 non-null float64
dtypes: float64(2), object(5)
memory usage: 7.7+ MB

2、数据切片分析

我们先来看一下variable,variable_full这两列的信息，

#数据切片分析
#看一下variable,variable_full这两列的信息，去掉重复项
data[['variable','variable_full']].drop_duplicates()#drop_duplicate方法是对DataFrame格式的数据，去除特定列下面的重复行。返回DataFrame格式的数据。
#total_area 国土面积(1000公顷)
#arable_land  可耕作面积
#permanent_crop_area 多年生作物面积
#cultivated_area  耕地面积
#percent_cultivated 耕地面积占比
#total_pop 总人口
#rural_pop 农村人口
#urban_pop 城市人口
#gdp 国内生产总值
#gdp_per_capita 人均国内生产总值
#agg_to_gdp 农业，增加国内生产总值
#human_dev_index  人类发展指数
#gender_inequal_index 性别不平等指数
#percent_undernourished 营养不良患病率
#avg_annual_rain_depth 长期平均年降水量
#national_rainfall_index  全国降雨指数

看一下统计了多少国家，

data.country.nunique()
输出结果
199

看一下有多少个时间周期，

data.time_period.nunique()#时间周期
输出结果
12

看一下时间周期有哪些，

time_period = data.time_period.unique()
time_period

我们看一下某一列某个指标的缺失值的个数，比如variable是total_area时缺失值的个数,

data[data.variable=='total_area'].value.isnull().sum()#缺失值个数
输出结果
220

我们通过几个维度来进行数据的分析：

• 横截面：一个时期内所有国家
• 时间序列：一个国家随着时间的推移
• 面板数据：所有国家随着时间的推移(作为数据给出)
• 地理空间：所有地理上相互联系的国家

我们按照上面的处理继续，现在我们想统计一下对于一个时间周期来说，不同国家在这个周期内的变化情况，

def time_slice(df,time_period):
    df = df[df.time_period == time_period]
    df = df.pivot(index = 'country',columns = 'variable',values='value')
    # 根据列对数据表进行重塑
    # index是重塑的新表的索引名称是什么
    # 第二个columns是重塑的新表的列名称是什么
    # values就是生成新列的值应该是多少
    df.columns.name=time_period
    return df

time_slice(data,time_period[0])

我们也可以按照国家分类，查看某个国家在不同时期的变化，

#按照国家分类，国家在不同时间周期内的变化
countries = data.country.unique()
def country_slice(df,country):
    df = df[df.country == country]
    df = df.pivot(index = 'variable',columns='time_period',values='value')
    df.index.name=country
    return df

country_slice(data,countries[40]).head(10)

我们还可以根据属性，查看不同国家在不同周期内的变化情况，

#对某一个属性或者变量进行分析
def variable_slice(df,variable):
    df = df[df.variable == variable]
    df = df.pivot(index = 'country',columns='time_period',values='value')
    df.index.name='country'
    return df

variable_slice(data,'total_pop').head(10)

我们还可以给定国家和指标，查看这个国家在这个指标上的变化情况，

#给定国家和指标，查看这个国家在这个指标上的变化情况
def time_series(df,country,variable):
    series = df[(df.country == country) & (df.variable == variable)]
    series = series.dropna()[['year_measured','value']]
    
    series.year_measured = series.year_measured.astype(int)#可用于转化dateframe某一列的数据类型
    series.set_index('year_measured',inplace=True)
    series.columns = [variable]
    return series
time_series(data,'Belarus','total_pop')

我们还有region(区域)没有查看，我们来看一下：

data.region.unique()

通过上图可以看出，区域太多，不便于观察，我们可以将一些区域进行合并。减少区域数量有助于模型评估，可以创建一个字典来查找新的，更简单的区域(亚洲，北美洲，南美洲，大洋洲)

simple_regions = {
    'World | Asia':'Asia',
       'Americas | Central America and Caribbean | Central America':'North America',
       'Americas | Central America and Caribbean | Greater Antilles':'North America',
       'Americas | Central America and Caribbean | Lesser Antilles and Bahamas':'North America',
       'Americas | Northern America | Northern America':'North America',
       'Americas | Northern America | Mexico':'North America',
       'Americas | Southern America | Guyana':'South America',
       'Americas | Southern America | Andean':'South America',
       'Americas | Southern America | Brazil':'South America',
       'Americas | Southern America | Southern America':'South America',
       'World | Africa':'Africa',
       'World | Europe':'Europe', 
       'World | Oceania':'Oceania'
}

data.region = data.region.apply(lambda x:simple_regions[x])

print(data.region.unique())
#输出结果
['Asia' 'North America' 'South America' 'Africa' 'Europe' 'Oceania']

我们来看一下数据变化，

data

#提取单个区域
def subregion(data,region):
    data = data[data.region == region]
    return data

subregion(data,'Asia').head(10)

3、单变量分析

紧接着上面的数据处理，我们重新导入一下包，这次有一些新包，

%matplotlib inline

#plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_context("poster",font_scale=1.3)
import folium  #画世界地图

import os,sys
import numpy as np
import pandas as pd

import pivottablejs
import missingno as msno#缺失值可视化
import pandas_profiling#可以以网页的形式展现给你数据总体概况

msno.matrix(recent,labels=True)#缺失值可视化

我们看一下水资源的情况，

#水资源总量
msno.matrix(variable_slice(data,'exploitable_total'),inline=False,sort='descending')
plt.xlabel('Time period')
plt.ylabel('Country')
plt.title('Missing total exploitable water resources data across countries and time  period \n \n \n')

通过上图可以看出只有一小部分国家报告了可利用的水资源总量，这些国家中只有极少数国家拥有最近一段时间的数据，我们将删除变量，因为这么少的数据点会导致很多问题。

data.loc[~data.variable.str.contains('exploitable'),:]

接下来我们看一下全国降雨指数，

# 全国降水指数（NRI）(毫米/每年)
msno.matrix(variable_slice(data,'national_rainfall_index'),inline=False,sort='descending')
plt.xlabel('Time period')
plt.ylabel('Country')
plt.title('Missing national rainfall index data across countries and time periods \n \n \n')

全国降雨在2002年以后不再报到，所以我们也删除这个数据，

data = data.loc[~(data.variable=='national_rainfall_index')]

我们单独拿出一个洲来进行分析，举例南美洲，我们来看一下数据的完整性，

#看美洲
# 先把数据筛选出来
north_america = subregion(data,'North America')

msno.matrix(msno.nullity_sort(time_slice(north_america,'2013-2017'),sort='descending').T,inline=False)#数据的完整性

msno.nullity_filter(country_slice(data,'Bahamas').T,filter='bottom',p=0.1)

#颜色越深，缺失值越严重
geo = r'../../数据集/农粮数据/world.json'
null_data = recent['agg_to_gdp'].notnull()*1
map = folium.Map(location=[48,-102],zoom_start=2)
map.choropleth(geo_data=geo,
              data=null_data,
              columns=['country','agg_to_gdp'],
              key_on='feature.properties.name',reset=True,
              fill_color='GnBu',fill_opacity=1,line_opacity=0.2,
              legend_name='Missing agricultural contribution to GDP data 2013-2017')

map

我们也可以指定不同的指标，

def plot_null_map(df,time_period,variable,legend_name=None):
    ts=time_slice(df,time_period).reset_index().copy()#不指明，从0开始
    ts[variable]=ts[variable].notnull()*1
    map = folium.Map(location=[48,-102],zoom_start=2)
    map.choropleth(geo_data=geo,
        data=ts,
        columns=['country',variable],
        key_on='feature.properties.name',reset=True,
        fill_color='GnBu',fill_opacity=1,line_opacity=0.2,
        legend_name=legend_name if legend_name else variable)
    return map

plot_null_map(data,'2013-2017','number_undernourished','Number undernourished id missing')

fig,ax = plt.subplots(figsize=(16,16))
sns.heatmap(data.groupby(['time_period','variable']).value.count().unstack().T,ax=ax)
plt.xticks(rotation=45)
plt.xlabel('Time period')
plt.ylabel('Variable')
plt.title('Number of countries with data reported for each variable over time')

接下来，我们使用pandas_profiling来对单变量以及多变量之间的关系进行统计一下，

pandas_profiling.ProfileReport(time_slice(data,'2013-2017'))

4、峰度与偏度

这里我们要计算的是，比如

• 均值，中位数，模式，四分位
• 标准差、方差、范围、间距范围
• 偏度、峰度
  首先我们先来看一下数据的描述，

recent[['total_pop','urban_pop','rural_pop']].describe().astype(int) #总人口，城镇人口，农村人口

recent.sort_values('rural_pop')[['total_pop','urban_pop','rural_pop']].head()

我们按照rural_pop从小到大进行排序，发现的确有几个国家的农村人口是负数，

time_series(data,'Qatar','total_pop').join(time_series(data,'Qatar','urban_pop')).join(time_series(data,'Qatar','rural_pop'))

人口数目是不可能小于0，所以这说明数据有问题，存在脏数据，如果做分析预测时，要注意将这些脏数据处理一下。

接下来我们看一下偏度，我们规定，

• 均值<中位数，定为左偏；
• 均值>中位数，定为有偏；

# 计算偏度
import scipy
recent[['total_pop','urban_pop','rural_pop']].apply(scipy.stats.skew)

正态分布的偏度应为零，负偏度表示左偏，正偏表示右偏。

偏度计算完后，我们计算一下峰度，峰度也是一个正态分布，峰度不能为负，只能是正数，越大说明越陡峭，

recent[['total_pop','urban_pop','rural_pop']].apply(scipy.stats.kurtosis)

接下来我们看一下，如果数据分布非常不均匀该怎么办呢，

fig,ax=plt.subplots(figsize=(12,6))
ax.hist(recent.total_pop.values,bins=50)
ax.set_xlabel('Total population')
ax.set_ylabel('Number of countries')
ax.set_title('Distribution of population of countries 2013-2017')

上图是2013-2017年国家总人数的分布，通过上图我们发现，人口量少于200000(不考虑单位)的国家非常多，人口大于1200000的国家非常少，如果我们需要建模的话，这种数据我们是不能要的。这个时候我们应该怎么办呢？

通常，遇到这种情况，使用log变换将其变为正常。对数变换是数据变换的一种常用方式，数据变换的目的在于使数据的呈现方式接近我们所希望的前提假设，从而更好的进行统计推断。

接下来，我们用log转换一下，并看一下它的偏度和峰值，

#偏度
recent[['total_pop']].apply(np.log).apply(scipy.stats.skew)

可以看出偏度下降了很多，减少了倾斜。

#峰度
recent[['total_pop']].apply(np.log).apply(scipy.stats.kurtosis)

可以发现峰度也下降了，接下来我们看一下经过log转换后的数据分布，

def plot_hist(df,variable,bins=20,xlabel=None,by=None,ylabel=None,title=None,logx=False,ax=None):
    if not ax:
        fig,ax=plt.subplots(figsize=(12,8))
    if logx:
        if df[variable].min()<=0:
            df[variable] = df[variable] -df[variable].min()+1
            print('Warning:data<=0 exists,data transformed by %0.2g before plotting' % (-df[variable].min()))
        bins = np.logspace(np.log10(df[variable].min()),np.log10(df[variable].max()),bins)
        ax.set_xscale("log")
    ax.hist(df[variable].dropna().values,bins=bins)
    
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title)
    
    return ax

plot_hist(recent,'total_pop',bins=25,logx=True,xlabel='Log of total population',ylabel='Number of countries',
         title='Distribution of total population of countries 2013-2017')

虽然数据还有一些偏度，但是明显好了很多，呈现的分布也比较标准。

5、数据的分析维度

首先我们先来看一下美国的人口总数随时间的变化，

plt.figure(figsize=(10,10))
plt.plot(time_series(data,'United States of America','total_pop'))
plt.xlabel('Year')
plt.ylabel('Population')
plt.title('United States population over time')

接下来，我们查看北美洲每个国家人口总数随着时间的变化，

plt.figure(figsize=(15, 10))
with sns.color_palette(sns.diverging_palette(220,280,s=85,l=25,n=23)):
    north_america = time_slice(subregion(data,'North America'),'1958-1962').sort_values('total_pop').index.tolist()
    for country in north_america:
        plt.plot(time_series(data,country,'total_pop'),label=country)
        plt.xlabel('Year')
        plt.ylabel('Population')
        plt.title('North American populations over time')

    plt.legend(loc=2,prop={'size':10})

这个时候我们发现，一些国家由于人口数量本身就少，所以整个图像显示的不明显，我们可以改变一下参照指标，那我们通过什么标准化？我们可以选择一个国家的最小、平均、中位数、最大值...或任何其他位置。那我们选择最小值，这样我们就能看到每个国家的起始人口上的增长。

plt.figure(figsize=(15, 10))
with sns.color_palette(sns.diverging_palette(220,280,s=85,l=25,n=23)):
    for country in north_america:
        ts=time_series(data,country,'total_pop')
        ts['norm_pop']=ts.total_pop/ts.total_pop.min()*100
        plt.plot(ts['norm_pop'],label=country)
        plt.xlabel('Year')
        plt.ylabel('Percent increase in population')
        plt.title('Percent increase in population from 1960 in North American countries')

    plt.legend(loc=2,prop={'size':10})

人口增长倍率图

我们也可以用热度图来展示，用颜色的深浅来比较大小关系，

north_america_pop = variable_slice(subregion(data,'North America'),'total_pop')
north_america_norm_pop = north_america_pop.div(north_america_pop.min(axis=1),axis=0)*100
north_america_norm_pop = north_america_norm_pop.loc[north_america]

fig,ax = plt.subplots(figsize=(16,12))
sns.heatmap(north_america_norm_pop,ax=ax,cmap=sns.light_palette((214,90,60),input="husl",as_cmap=True))
plt.xticks(rotation=45)
plt.xlabel('Time period')
plt.ylabel('Country,ordered by population in 1960(<-greatest to least->)')
plt.title('Percent increase in population from 1960')

接下来我们分析一下水资源的分布情况，

plot_hist(recent,'total_renewable',bins=50,
         xlabel='Total renewable water resources (10^9 m^3/yr$)',
         ylabel='Number of countries',
         title='Distribution of total renewable water resources,2013-2017')

我们可以进行一下log转换，

plot_hist(recent,'total_renewable',bins=50,
         xlabel='Total renewable water resources (10^9 m^3/yr$)',
         ylabel='Number of countries',logx=True,
         title='Distribution of total renewable water resources,2013-2017')

我们用热度图画一下，

north_america_renew = variable_slice(subregion(data,'North America'),'total_renewable')

fig,ax = plt.subplots(figsize=(16,12))
sns.heatmap(north_america_renew,ax=ax,cmap=sns.light_palette((214,90,60),input="husl",as_cmap=True))
plt.xticks(rotation=45)
plt.xlabel('Time period')
plt.ylabel('Country,ordered by Total renewable water resources in 1960(<-greatest to least->)')
plt.title('Total renewable water resources increase in population from 1960')

6、变量关系可视化展示

连续值可以画成散点图，方便观看，
我们来看一下随着季节变化，人均GDP的变化情况，

data=data.loc[~data.variable.str.contains('exploitable')]
data=data.loc[~(data.variable=='national_rainfall_index')]

plt.figure(figsize=(8,8))
plt.scatter(recent.seasonal_variability,recent.gdp_per_capita)
plt.xlabel('Seasonal variability')
plt.ylabel('GDP per capita ($USD/person)')

def plot_scatter(df,x,y,bins=20,xlabel=None,ylabel=None,title=None,ax=None,logx=False,logy=False):
    if not ax:
        fig,ax=plt.subplots(figsize=(12,8))
    colors = mpl.rcParams['axes.prop_cycle'].by_key()['color']
    if by:
        groups = df.groupby(by)
        for j,(name,group) in enumerate(groups):
            ax.scatter(group[x],group[y],color=colors[j],label=name)
        ax.legend()
    else:
        ax.scatter(df[x],df[y],color=colors[0])
    if logx:
        ax.set_xscale('log')
    if logy:
        ax.set_yscale('log')
    ax.set_xlabel(xlabel if xlabel else x)
    ax.set_ylabel(ylabel if ylabel else y)
    if title:
        ax.set_title(title)
    
    return ax

svr = [recent.seasonal_variability.min(),recent.seasonal_variability.max()]
gdpr = [recent.gdp_per_capita.min(),recent.gdp_per_capita.max()]
gdpbins = np.logspace(*np.log10(gdpr),25)

g = sns.JointGrid(x="seasonal_variability",y="gdp_per_capita",data=recent,ylim=gdpr)
g.ax_marg_x.hist(recent.seasonal_variability,range=svr)
g.ax_marg_y.hist(recent.gdp_per_capita,range=gdpr,bins=gdpbins,orientation="horizontal")
g.plot_joint(plt.hexbin,gridsize=25)
ax = g.ax_joint

g.fig.set_figheight(8)
g.fig.set_figwidth(9)

相关程度：
相关度量两个变量之间的线性关系的强度，我们可以用相关性来识别变量。

现在我们单独拿出来一个指标分析是什么因素与人均GDP的变化有关系，正相关就是积极影响，负相关就是消极影响。

recent_corr = recent.corr().loc['gdp_per_capita'].drop(['gdp','gdp_per_capita'])

def conditional_bar(series,bar_colors=None,color_labels=None,figsize=(13,24),
                   xlabel=None,by=None,ylabel=None,title=None):
    fig,ax = plt.subplots(figsize=figsize)
    if not bar_colors:
        bar_colors = mpl.rcParams['axes.prop_cycle'].by_key()['color'][0]
    plt.barh(range(len(series)),series.values,color=bar_colors)
    plt.xlabel('' if not xlabel else xlabel)
    plt.ylabel('' if not ylabel else ylabel)
    plt.yticks(range(len(series)),series.index.tolist())
    plt.title('' if not title else title)
    plt.ylim([-1,len(series)])
    if color_labels:
        for col,lab in color_labels.items():
            plt.plot([],linestyle='',marker='s',c=col,label=lab)
        lines,labels=ax.get_legend_handles_labels()
        ax.legend(lines[-len(color_labels.keys()):],labels[-len(color_labels.keys()):],loc='upper right')
    plt.close()
    return fig

bar_colors = ['#0055A7' if x else '#2C3E4F' for x in list(recent_corr.values<0)]
color_labels = {'#0055A7':'Negative correlation','#2C3E4F':'Positive correlation'}

conditional_bar(recent_corr.apply(np.abs),bar_colors,color_labels,
               title='Magnitude of correlation with GDP per capita,2013-2017',
               xlabel='|Correlation|')

当我们在画图的时候也可以考虑一下利用bined设置一下区间，比如说连续值我们可以分成几个区间进行分析，这里我们以人均GDP的数量来进行分析，我们可以将人均GDP的数据映射到不同的区间，比如人均GDP比较低，比较落后的国家，以及人均GDP比较高，比较发达的国家，这个也是我们经常需要的操作，

plot_hist(recent,'gdp_per_capita',xlabel='GDP per capita($)',
         ylabel='Number of countries',
         title='Distribution of GDP per capita,2013-2017')

做一下log变换，这里是25个bin

plot_hist(recent,'gdp_per_capita',xlabel='GDP per capita($)',logx=True,
         ylabel='Number of countries',bins=25,
         title='Distribution of GDP per capita,2013-2017')

我们指定一下分割的标准，

capita_bins = ['Very low','Low','Medium','High','Very high']
recent['gdp_bin'] = pd.qcut(recent.gdp_per_capita,5,capita_bins)
bin_ranges = pd.qcut(recent.gdp_per_capita,5).unique()

def plot_hist_div(df,variable,bins=None,xlabel=None,by=None,ylabel=None,title=None,logx=False,ax=None):
    if not ax:
        fig,ax=plt.subplots(figsize=(12,8))
    if logx:
        bins = np.logspace(np.log10(df[variable].min()),np.log10(df[variable].max()),bins)#返回数以对数刻度均匀分布。
        ax.set_xscale("log")
    if by:
        if type(df[by].unique()) == pd.Categorical:
            cats = df[by].unique().categories.tolist()
        else:
            cats = df[by].unique().tolist()
        for cat in cats:
            to_plot = df[df[by] == cat][variable].dropna()
            ax.hist(to_plot,bins=bins)
    else:
        ax.hist(df[variable].dropna().values,bins=bins)
    
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title)
    
    return ax

plot_hist_div(recent,'gdp_per_capita',xlabel='GDP per capita($)',logx=True,
         ylabel='Number of countries',bins=25,by='gdp_bin',
         title='Distribution of log GDP per capita,2013-2017')

我们还可以看一下人均GDP较低，落后国家的内部数据，下面我们看一下内部数据分布情况，用boxplot进行画图，

recent[['gdp_bin','total_pop_access_drinking']].boxplot(by='gdp_bin',figsize=(10,10))
plt.title('Distribution of percent of total population with access to drinking water across gdp per capita categories')
plt.xlabel('GDP per capita quintile')
plt.ylabel('Total population of country')

对于这部分的分布，我们还可以统计看一下其他指标，如下图所示，我们还可以看一下洪水的统计信息，

def mult_boxplots(df,variable,category,
                 xlabel=None,ylabel=None,
                 title=None,ylim=None):
    df[[variable,category]].boxplot(by=category,figsize=(10,10))
    
    if xlabel:
        plt.xlabel(xlabel)
    if ylabel:
        plt.ylabel(ylabel)
    if title:
        plt.title(title)
    if ylim:
        plt.ylim(ylim)

mult_boxplots(recent,'flood_occurence','gdp_bin',
             xlabel='GDP per capita quintile')

我们现在总结一下上面用的函数，为了以后方便使用，我们将用到的函数在统计一下，

#某一个时间区域内，各个国家与各个指标之间的关系
def time_slice(df,time_period):
    df = df[df.time_period == time_period]
    df = df.pivot(index = 'country',columns = 'variable',values='value')
    # 根据列对数据表进行重塑
    # index是重塑的新表的索引名称是什么
    # 第二个columns是重塑的新表的列名称是什么
    # values就是生成新列的值应该是多少
    df.columns.name=time_period
    return df

#指定国家，查看各个指标随时间的变化
countries = data.country.unique()
def country_slice(df,country):
    df = df[df.country == country]
    df = df.pivot(index = 'variable',columns='time_period',values='value')
    df.index.name=country
    return df


#指定指标，查看各个国家随着时间的变化
def variable_slice(df,variable):
    df = df[df.variable == variable]
    df = df.pivot(index = 'country',columns='time_period',values='value')
    df.index.name='country'
    return df


#给定国家和指标，查看这个国家在这个指标上的变化情况
def time_series(df,country,variable):
    series = df[(df.country == country) & (df.variable == variable)]
    series = series.dropna()[['year_measured','value']]
    
    series.year_measured = series.year_measured.astype(int)#可用于转化dateframe某一列的数据类型
    series.set_index('year_measured',inplace=True)
    series.columns = [variable]
    return series

#将国家按照各大洲分类
simple_regions = {
    'World | Asia':'Asia',
       'Americas | Central America and Caribbean | Central America':'North America',
       'Americas | Central America and Caribbean | Greater Antilles':'North America',
       'Americas | Central America and Caribbean | Lesser Antilles and Bahamas':'North America',
       'Americas | Northern America | Northern America':'North America',
       'Americas | Northern America | Mexico':'North America',
       'Americas | Southern America | Guyana':'South America',
       'Americas | Southern America | Andean':'South America',
       'Americas | Southern America | Brazil':'South America',
       'Americas | Southern America | Southern America':'South America',
       'World | Africa':'Africa',
       'World | Europe':'Europe', 
       'World | Oceania':'Oceania'
}
data.region = data.region.apply(lambda x:simple_regions[x])


#提取单个区域
def subregion(data,region):
    data = data[data.region == region]
    return data

#查看各个国家总人口的分布情况(数据经过log转换)
def plot_hist(df,variable,bins=20,xlabel=None,by=None,ylabel=None,title=None,logx=False,ax=None):
    if not ax:
        fig,ax=plt.subplots(figsize=(12,8))
    if logx:
        if df[variable].min()<=0:
            df[variable] = df[variable] -df[variable].min()+1
            print('Warning:data<=0 exists,data transformed by %0.2g before plotting' % (-df[variable].min()))
        bins = np.logspace(np.log10(df[variable].min()),np.log10(df[variable].max()),bins)#返回数以对数刻度均匀分布。
        ax.set_xscale("log")
    ax.hist(df[variable].dropna().values,bins=bins)
    
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title)
    
    return ax


#提取某个洲的各个国家，查看总人口随时间的变化
plt.figure(figsize=(15, 10))
with sns.color_palette(sns.diverging_palette(220,280,s=85,l=25,n=23)):
    north_america = time_slice(subregion(data,'North America'),'1958-1962').sort_values('total_pop').index.tolist()
    for country in north_america:
        plt.plot(time_series(data,country,'total_pop'),label=country)
        plt.xlabel('Year')
        plt.ylabel('Population')
        plt.title('North American populations over time')

    plt.legend(loc=2,prop={'size':10})


#指定人口基数，查看各个国家在这个基数上人口变化的倍率
plt.figure(figsize=(15, 10))
with sns.color_palette(sns.diverging_palette(220,280,s=85,l=25,n=23)):
    for country in north_america:
        ts=time_series(data,country,'total_pop')
        ts['norm_pop']=ts.total_pop/ts.total_pop.min()*100
        plt.plot(ts['norm_pop'],label=country)
        plt.xlabel('Year')
        plt.ylabel('Percent increase in population')
        plt.title('Percent increase in population from 1960 in North American countries')

    plt.legend(loc=2,prop={'size':10})


#相关性的条形图
def conditional_bar(series,bar_colors=None,color_labels=None,figsize=(13,24),
                   xlabel=None,by=None,ylabel=None,title=None):
    fig,ax = plt.subplots(figsize=figsize)
    if not bar_colors:
        bar_colors = mpl.rcParams['axes.prop_cycle'].by_key()['color'][0]
    plt.barh(range(len(series)),series.values,color=bar_colors)
    plt.xlabel('' if not xlabel else xlabel)
    plt.ylabel('' if not ylabel else ylabel)
    plt.yticks(range(len(series)),series.index.tolist())
    plt.title('' if not title else title)
    plt.ylim([-1,len(series)])
    if color_labels:
        for col,lab in color_labels.items():
            plt.plot([],linestyle='',marker='s',c=col,label=lab)
        lines,labels=ax.get_legend_handles_labels()
        ax.legend(lines[-len(color_labels.keys()):],labels[-len(color_labels.keys()):],loc='upper right')
    plt.close()
    return fig


#将数据分开不同的区域进行相关性分析
def plot_hist_div(df,variable,bins=None,xlabel=None,by=None,ylabel=None,title=None,logx=False,ax=None):
    if not ax:
        fig,ax=plt.subplots(figsize=(12,8))
    if logx:
        bins = np.logspace(np.log10(df[variable].min()),np.log10(df[variable].max()),bins)#返回数以对数刻度均匀分布。
        ax.set_xscale("log")
    if by:
        if type(df[by].unique()) == pd.Categorical:
            cats = df[by].unique().categories.tolist()
        else:
            cats = df[by].unique().tolist()
        for cat in cats:
            to_plot = df[df[by] == cat][variable].dropna()
            ax.hist(to_plot,bins=bins)
    else:
        ax.hist(df[variable].dropna().values,bins=bins)
    
    if xlabel:
        ax.set_xlabel(xlabel)
    if ylabel:
        ax.set_ylabel(ylabel)
    if title:
        ax.set_title(title)
    
    return ax


#指标的不同类别进行boxplot分析
def mult_boxplots(df,variable,category,
                 xlabel=None,ylabel=None,
                 title=None,ylim=None):
    df[[variable,category]].boxplot(by=category,figsize=(10,10))
    
    if xlabel:
        plt.xlabel(xlabel)
    if ylabel:
        plt.ylabel(ylabel)
    if title:
        plt.title(title)
    if ylim:
        plt.ylim(ylim)

关于农粮组织数据的探索性分析到这里就结束了，希望大家能动手做一下。