利用Python可视化全球二氧化碳排放

大家好，近几年碳的问题一直很火，那么关于碳数据的收集就比较重要，刚好前不久见到一个全球实时碳数据的网站，里面包含碳的数据，网站包含一些主要国家碳排放数据。网站链接：全球实时碳数据 Carbon Monitor

关于我们如何绘制地图，刚好有个国外大神Adam Symington画了很多好看的地图，其中就有二氧化碳排放的地图。我们借鉴他的思路，选择自己想要表达的效果。接下来，我们开始对数据进行可视化。

一.读取二氧化碳排放数据

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# author：GuaGua
# datetime：2022/5/18 11:15
# env:python3.8
import pandas as pd
#读取数据
data= pd.read_csv("v50_CO2_excl_short-cycle_org_C_2018.txt", delimiter=';')
data

	lat	lon	emission
0	88.1	-50.7	11.0301
1	88.1	-50.6	15.7066
2	88.1	-50.5	19.1943
3	88.1	-50.4	22.0621
4	88.1	-50.3	24.5536
...	...	...	...
2191880	-77.8	166.5	11459.2000
2191881	-77.8	166.6	694.4980
2191882	-77.9	166.3	2133.1000
2191883	-77.9	166.4	1389.0000
2191884	-77.9	166.5	7887.5100

2191885 rows × 3 columns

二.可视化画图

#首先进行简单可视化
import matplotlib.pyplot as plt
#设置字体
plt.rc('font',family='Times New Roman') 
fig = plt.figure(figsize=(10,5))
pic=plt.scatter(data['lon'], data['lat'],c=data['emission'], s=0.05,  edgecolors='none')
plt.colorbar(pic)
plt.show()

三.对图形进行优化

#对数颜色映射的一个小例子，帮助大家理解下面颜色映射，为后续做好准备
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
#首先构造一个指数函数
x=np.linspace(0,10,100)
y=np.linspace(0,10,100)
matX=np.matrix(np.exp(2*x))
matY=np.matrix(np.exp(3*y)).T
Z=np.array(matY @ matX)
fig, ax = plt.subplots(3, 1,figsize=(10,8))
pcm=ax[0].contourf(x,y,Z,cmap='jet',levels=100)
fig.colorbar(pcm,ax=ax[0])
#matplotlib.colors.LogNorm(vmin=None, vmax=None, clip=False)
#将给定值对应于对数化的0-1的范围，其中vmin和vmax定义的是用来scale的值的区间，如果这两个值没有给，那么默认是最大值和最小值。
#我们可以看见绝大部分都是蓝色，只有尖角的一小部分有红色，默认的颜色分割方式把它等距的分割开，这个情况下，我们使用对数尺度的颜色映射
pcm=ax[1].contourf(x,y,Z,cmap='jet',norm=colors.LogNorm())
fig.colorbar(pcm,ax=ax[1])
logmax=np.log10(Z.max())
logmin=np.log10(Z.min())
#按照对数均匀生成100个等级。三种颜色映射对比
lvls=np.logspace(logmin,logmax,100)
pcm=ax[2].contourf(x,y,Z,cmap='jet',levels=lvls,norm=colors.LogNorm())
fig.colorbar(pcm,ax=ax[2])
plt.show()

#原始颜色对数映射。对比更加明显
from matplotlib import colors
fig = plt.figure(figsize=(10,5))
pic = plt.scatter(data['lon'], data['lat'],c=data['emission'],s=0.05, edgecolors='none',  norm=colors.LogNorm())
plt.colorbar(pic)
plt.show()

由于matplotlib提供的颜色映射表是有限的，所以我们还需要借助外部的库包提供额外的颜色映射表。大气科学与海洋科学常用的外部颜色库包常为cmaps，cmaps库包是将NCL平台的颜色移植到python。

import cmaps
import numpy as np
import inspect
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rc('text', usetex=False)
#展示所有的颜色映射
def list_cmaps():
    attributes = inspect.getmembers(cmaps, lambda _: not (inspect.isroutine(_)))
    colors = [_[0] for _ in attributes if
              not (_[0].startswith('__') and _[0].endswith('__'))]
    return colors


if __name__ == '__main__':
    color = list_cmaps()

    a = np.outer(np.arange(0, 1, 0.001), np.ones(10))
    plt.figure(figsize=(20, 20))
    plt.subplots_adjust(top=0.95, bottom=0.05, left=0.01, right=0.99)
    ncmaps = len(color)
    nrows = 8
    for i, k in enumerate(color):
        plt.subplot(nrows, ncmaps // nrows + 1, i + 1)
        plt.axis('off')
        plt.imshow(a, aspect='auto', cmap=getattr(cmaps, k), origin='lower')
        plt.title(k, rotation=90, fontsize=10)
        plt.title(k, fontsize=10)
    plt.savefig('colormaps.png', dpi=300)

import cmaps   
#调整颜色,根据上面的色带选择自己喜欢的颜色
fig, ax = plt.subplots(2, 1,figsize=(10,10))
#在使用cmaps时，只需引入库包，然后填写颜色名称即可,选用两种色带进行对比看看
cmap1=cmaps.BlueWhiteOrangeRed
cmap2=cmaps.MPL_RdYlGn_r
pic = ax[0].scatter(data['lon'], data['lat'],c=data['emission'],s=0.05, edgecolors='none',  norm=colors.LogNorm(),cmap=cmap1)
fig.colorbar(pic,ax=ax[0])
pic = ax[1].scatter(data['lon'], data['lat'],c=data['emission'],s=0.05, edgecolors='none',  norm=colors.LogNorm(),cmap=cmap2)
fig.colorbar(pic,ax=ax[1])
plt.show()

四.对图形进行投影

import geopandas as gpd
from shapely.geometry import Point
geometry = [Point(xy) for xy in zip(data['lon'], data['lat'])]# # 需要修改为对应的经纬度字段
geodata = gpd.GeoDataFrame(data, crs="EPSG:4326", geometry=geometry)# 指定坐标系
geodata

import matplotlib.pyplot as plt                                          
import cartopy.crs as ccrs                                               
from matplotlib import cm
from matplotlib.colors import ListedColormap
import matplotlib as mpl   
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
plt.rc('font',family='Times New Roman') 
fig=plt.figure(figsize=(6,4),dpi=200)
#ax首先得定义坐标系
ax=fig.add_axes([0,0,1,1],projection=ccrs.PlateCarree())
ax = geodata.plot(ax=ax, column='emission', transform=ccrs.PlateCarree(),
                  cmap=cmap1, norm=colors.LogNorm(), s=0.05, edgecolors='none')
ax.set_title('2018 CO$_2$ Emissions',fontsize=15,fontweight ='bold')
#添加文本信息
text = ax.text(0.0, 0.02, "By GuaGua",
              size=10,
              fontweight ='bold',
              transform = ax.transAxes)
#gridlines设置经纬度标签
gl=ax.gridlines(draw_labels=True,linestyle=":",linewidth=0.3,color='red')
gl.top_labels=False #关闭上部经纬标签                                  
gl.right_labels=False#关闭右边经纬标签  
gl.xformatter = LONGITUDE_FORMATTER  #使横坐标转化为经纬度格式            
gl.yformatter = LATITUDE_FORMATTER                                        
gl.xlocator=mticker.FixedLocator(np.arange(-180,180,20))      
gl.ylocator=mticker.FixedLocator(np.arange(-90,90,20)) 
#修改经纬度字体大小 
gl.xlabel_style={'size':8, 'weight': 'bold'}                            
gl.ylabel_style={'size':8, 'weight': 'bold'}

#使用Robinson投影，默认中央子午线为中心
import cartopy.crs as ccrs 
import matplotlib.ticker as mticker
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10,5), subplot_kw={'projection': ccrs.Robinson()})
ax = geodata.plot(ax=ax, column='emission', transform=ccrs.PlateCarree(),
                  cmap=cmap2, norm=colors.LogNorm(), s=0.05, edgecolors='none')
#添加文本信息
text = ax.text(0.0, 0.02, "By GuaGua",
              size=10,
              fontweight ='bold',
              transform = ax.transAxes)
#添加格网
g2=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g2.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
g2.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
#设置经纬度网格的间隔
g2.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g2.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改经纬度字体大小 
g2.xlabel_style={'size':10,'weight': 'bold'}                            
g2.ylabel_style={'size':10,'weight': 'bold'}
#添加标题
ax.set_title('2018 CO$_2$ Emissions',
              fontsize=15,
              pad=20,
              fontweight ='bold')
plt.show()

五.色彩映射

#发现前面颜色对比不够明显，我们重新设置colorbar，截取部分colormap
from matplotlib import cm
from matplotlib.colors import ListedColormap
import numpy as np
import matplotlib as mpl   
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(10,8),dpi=300)
ax1=fig.add_axes([0,0,1,0.05])                                                 
ax2=fig.add_axes([0,0.1,1,0.05])  
ax3=fig.add_axes([0,0.2,1,0.05])  
#第一个 
our_cmap = cm.get_cmap('afmhot', 11)
newcolors = our_cmap(np.linspace(0, 1, 11))#分片操作，生成0到1的11个数据间隔的数组   
our_cmap = ListedColormap(newcolors[::1])#重构为新的colormap
#第二个                                    
newcolors_1=cmap2(np.linspace(0,1,11))#分片操作，生成0到1的11个数据间隔的数组                                  
newcmap=ListedColormap(newcolors_1[::1])  #重构为新的colormap
#第三个
newcolors_2=cmap1(np.linspace(0,1,11))#分片操作，生成0到1的11个数据间隔的数组                                  
newcmap2=ListedColormap(newcolors_2[::1])  #重构为新的colormap
###########################################
bounds = [0.0, 0.06, 6, 60, 600, 3000, 6000, 24000, 45000, 120000]
norm = colors.BoundaryNorm(bounds, our_cmap.N)#基于离散间隔创建颜色图
norm1 = colors.BoundaryNorm(bounds, newcmap.N)#基于离散间隔创建颜色图
norm2 = colors.BoundaryNorm(bounds, newcmap2.N)#基于离散间隔创建颜色图
fc1=fig.colorbar(mpl.cm.ScalarMappable(norm=norm,cmap=our_cmap),              
                cax=ax1,                                                      
                orientation='horizontal',   
                label='our cmap',
                ticks=bounds,                                                    
                extend='both')
fc2=fig.colorbar(mpl.cm.ScalarMappable(norm=norm1,cmap=newcmap),              
                cax=ax2,                                                      
                orientation='horizontal',   
                label='newcmap',
                ticks=bounds,                                                    
                extend='both')
fc3=fig.colorbar(mpl.cm.ScalarMappable(norm=norm2,cmap=newcmap2),              
                cax=ax3,                                                      
                orientation='horizontal',   
                label='newcmap2',
                ticks=bounds,                                                    
                extend='both')

fig, ax = plt.subplots(figsize=(10,5), subplot_kw={'projection': ccrs.Robinson()})
ax = geodata.plot(ax=ax, column='emission', transform=ccrs.PlateCarree(),
                  cmap=newcmap, norm=norm1, s=0.05, alpha=1, edgecolors='none')
#添加文本信息
text = ax.text(0.0, 0.02, "By GuaGua",
              size=10,
              fontweight ='bold',
              transform = ax.transAxes)
#添加格网
g3=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g3.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
g3.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
#设置经纬度网格的间隔
g3.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g3.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改经纬度字体大小 
g3.xlabel_style={'size':10,'weight': 'bold'}                            
g3.ylabel_style={'size':10,'weight': 'bold'}
#添加标题
ax.set_title('2018 CO$_2$ Emissions',
              fontsize=15,
              fontweight ='bold',
              pad=20)
plt.savefig('co2_emissions_noclorbar.png', dpi=300,bbox_inches='tight')
plt.show()

六.添加图例

fig, ax = plt.subplots(figsize=(10,5),subplot_kw={'projection': ccrs.Robinson()})
ax = geodata.plot(ax=ax, column='emission', transform=ccrs.PlateCarree(),
                  cmap=newcmap, norm=norm1, s=0.05, alpha=1, edgecolors='none')
#添加文本信息
text= ax.text(0.0, 0.02, "By GuaGua",
              size=10,
              fontweight ='bold',
              transform = ax.transAxes)
#添加格网
g3=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g3.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
g3.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
#设置经纬度网格的间隔
g3.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g3.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改经纬度字体大小 
g3.xlabel_style={'size':10,'weight': 'bold'}                            
g3.ylabel_style={'size':10,'weight': 'bold'}
#添加标题
ax.set_title('2018 CO$_2$ Emissions',
              fontsize=15,
              fontweight ='bold',
              pad=20)
#设置图例
fig = ax.get_figure()
cax = fig.add_axes([0.36, 0.16, 0.33, 0.01])
sm = plt.cm.ScalarMappable(cmap=newcmap, norm=norm1)
cb = fig.colorbar(sm, cax=cax, orientation="horizontal", pad=0.2, format='%.1e',
                  ticks=[0.03, 3, 33, 330, 1800, 4500, 15000, 34500, 82500],
                  drawedges=True)
cb.outline.set_visible(False)
cb.ax.tick_params(labelsize=6, width=0.5, length=0.5)
cbytick_obj = plt.getp(cb.ax, 'xticklabels' ) 
plt.setp(cbytick_obj)
cb.ax.set_xlabel('CO$_2$ Tons/Year', fontsize=6,  labelpad=-20)
# plt.savefig('co2_emissions.jpg', dpi=300, bbox_inches='tight')
plt.show()