編程分享|利用Python可視化全球二氧化碳排放

原創(chuàng) GuaGua 農(nóng)業(yè)與土地遙感 2022-05-19 08:30 發(fā)表于美國(guó)

背景

近幾年碳的問(wèn)題一直很火，那么關(guān)于碳數(shù)據(jù)的收集就比較重要，剛好前不久見(jiàn)到一個(gè)全球?qū)崟r(shí)碳數(shù)據(jù)的網(wǎng)站，里面包含碳的數(shù)據(jù)，網(wǎng)站包含一些主要國(guó)家碳排放數(shù)據(jù)，碳排放數(shù)據(jù)涉及到的來(lái)源有： “電力”、“地面運(yùn)輸”、“工業(yè)”、“居民消費(fèi)”、“國(guó)內(nèi)航空”、“國(guó)際航空”。

網(wǎng)站鏈接：https:///

此網(wǎng)站不僅數(shù)據(jù)好，而且圖也很漂亮，剛打開(kāi)首頁(yè)，就被圖吸引住了。

那么我們?nèi)绾卫L制地圖呢？剛好有個(gè)國(guó)外大神Adam Symington畫(huà)了很多好看的地圖，其中就有二氧化碳排放的地圖。我們借鑒模仿大神的思路，選擇自己想要表達(dá)的效果。

Adam Symington作品

首先是獲取數(shù)據(jù)，數(shù)據(jù)來(lái)源同大神作品的數(shù)據(jù)，來(lái)源于另一個(gè)碳數(shù)據(jù)網(wǎng)站：https://edgar.jrc.ec./  此網(wǎng)站有許多數(shù)據(jù)集可供您選擇，不僅僅有碳排放數(shù)據(jù)。

在本文中，我們只研究 2018 年全球的二氧化碳排放量，其他數(shù)據(jù)不進(jìn)行研究，感興趣的可以去原網(wǎng)站探索。數(shù)據(jù)的描述和原理見(jiàn)原網(wǎng)站，本文不再做贅述，直接開(kāi)始對(duì)數(shù)據(jù)進(jìn)行可視化。

1.讀取二氧化碳排放數(shù)據(jù)

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

2191885 rows × 3 columns

2.可視化畫(huà)圖

#首先進(jìn)行簡(jiǎn)單可視化
import matplotlib.pyplot as plt
#設(shè)置字體
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()

3.對(duì)圖形進(jìn)行優(yōu)化一下

#對(duì)數(shù)顏色映射的一個(gè)小例子，幫助大家理解下面顏色映射，為后續(xù)做好準(zhǔn)備
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
#首先構(gòu)造一個(gè)指數(shù)函數(shù)
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)
#將給定值對(duì)應(yīng)于對(duì)數(shù)化的0-1的范圍，其中vmin和vmax定義的是用來(lái)scale的值的區(qū)間，如果這兩個(gè)值沒(méi)有給，那么默認(rèn)是最大值和最小值。
#我們可以看見(jiàn)絕大部分都是藍(lán)色，只有尖角的一小部分有紅色，默認(rèn)的顏色分割方式把它等距的分割開(kāi)，這個(gè)情況下，我們使用對(duì)數(shù)尺度的顏色映射
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())
#按照對(duì)數(shù)均勻生成100個(gè)等級(jí)。三種顏色映射對(duì)比
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()

#原始顏色對(duì)數(shù)映射。對(duì)比更加明顯
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提供的顏色映射表是有限的，所以我們還需要借助外部的庫(kù)包提供額外的顏色映射表。大氣科學(xué)與海洋科學(xué)常用的外部顏色庫(kù)包常為cmaps，cmaps庫(kù)包是將NCL平臺(tái)的顏色移植到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   
#調(diào)整顏色,根據(jù)上面的色帶選擇自己喜歡的顏色
fig, ax = plt.subplots(2, 1,figsize=(10,10))
#在使用cmaps時(shí)，只需引入庫(kù)包，然后填寫(xiě)顏色名稱(chēng)即可,選用兩種色帶進(jìn)行對(duì)比看看
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()

4.對(duì)圖形進(jìn)行投影

import geopandas as gpd
from shapely.geometry import Point
geometry = [Point(xy) for xy in zip(data['lon'], data['lat'])]# # 需要修改為對(duì)應(yīng)的經(jīng)緯度字段
geodata = gpd.GeoDataFrame(data, crs="EPSG:4326", geometry=geometry)# 指定坐標(biāo)系
geodata

2191885 rows × 4 columns

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首先得定義坐標(biāo)系
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設(shè)置經(jīng)緯度標(biāo)簽
gl=ax.gridlines(draw_labels=True,linestyle=":",linewidth=0.3,color='red')
gl.top_labels=False #關(guān)閉上部經(jīng)緯標(biāo)簽                                  
gl.right_labels=False#關(guān)閉右邊經(jīng)緯標(biāo)簽  
gl.xformatter = LONGITUDE_FORMATTER  #使橫坐標(biāo)轉(zhuǎn)化為經(jīng)緯度格式            
gl.yformatter = LATITUDE_FORMATTER                                        
gl.xlocator=mticker.FixedLocator(np.arange(-180,180,20))      
gl.ylocator=mticker.FixedLocator(np.arange(-90,90,20)) 
#修改經(jīng)緯度字體大小 
gl.xlabel_style={'size':8, 'weight': 'bold'}                            
gl.ylabel_style={'size':8, 'weight': 'bold'}

#使用Robinson投影，默認(rèn)中央子午線(xiàn)為中心
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)
#添加格網(wǎng)
g2=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g2.xformatter = LONGITUDE_FORMATTER  # x軸設(shè)為經(jīng)度的格式
g2.yformatter = LATITUDE_FORMATTER  # y軸設(shè)為緯度的格式
#設(shè)置經(jīng)緯度網(wǎng)格的間隔
g2.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g2.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改經(jīng)緯度字體大小 
g2.xlabel_style={'size':10,'weight': 'bold'}                            
g2.ylabel_style={'size':10,'weight': 'bold'}
#添加標(biāo)題
ax.set_title('2018 CO$_2$ Emissions',
              fontsize=15,
              pad=20,
              fontweight ='bold')
plt.show()

#使用Robinson投影，以中國(guó)為中心
fig, ax = plt.subplots(figsize=(10,5),subplot_kw={'projection': ccrs.Robinson(central_longitude=130)})
ax = geodata.plot(ax=ax, column='emission', transform=ccrs.PlateCarree(),
                  cmap='afmhot', norm=colors.LogNorm(), s=0.05, edgecolors='none')
#添加文本信息
text = ax.text(0.0, 0.02, "By GuaGua",
              size=10,
              fontweight ='bold',
              transform = ax.transAxes)
#添加格網(wǎng)
g3=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g3.xformatter = LONGITUDE_FORMATTER  # x軸設(shè)為經(jīng)度的格式
g3.yformatter = LATITUDE_FORMATTER  # y軸設(shè)為緯度的格式
#設(shè)置經(jīng)緯度網(wǎng)格的間隔
g3.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g3.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改經(jīng)緯度字體大小 
g3.xlabel_style={'size':10,'weight': 'bold'}                            
g3.ylabel_style={'size':10,'weight': 'bold'}
#添加標(biāo)題
ax.set_title('2018 CO$_2$ Emissions',
              fontsize=15,
              fontweight ='bold',
              pad=20)
plt.show()

5.色彩映射

#發(fā)現(xiàn)前面顏色對(duì)比不夠明顯，我們重新設(shè)置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])  
#第一個(gè) 
our_cmap = cm.get_cmap('afmhot', 11)
newcolors = our_cmap(np.linspace(0, 1, 11))#分片操作，生成0到1的11個(gè)數(shù)據(jù)間隔的數(shù)組   
our_cmap = ListedColormap(newcolors[::1])#重構(gòu)為新的colormap
#第二個(gè)                                    
newcolors_1=cmap2(np.linspace(0,1,11))#分片操作，生成0到1的11個(gè)數(shù)據(jù)間隔的數(shù)組                                  
newcmap=ListedColormap(newcolors_1[::1])  #重構(gòu)為新的colormap
#第三個(gè)
newcolors_2=cmap1(np.linspace(0,1,11))#分片操作，生成0到1的11個(gè)數(shù)據(jù)間隔的數(shù)組                                  
newcmap2=ListedColormap(newcolors_2[::1])  #重構(gòu)為新的colormap
###########################################
bounds = [0.0, 0.06, 6, 60, 600, 3000, 6000, 24000, 45000, 120000]
norm = colors.BoundaryNorm(bounds, our_cmap.N)#基于離散間隔創(chuàng)建顏色圖
norm1 = colors.BoundaryNorm(bounds, newcmap.N)#基于離散間隔創(chuàng)建顏色圖
norm2 = colors.BoundaryNorm(bounds, newcmap2.N)#基于離散間隔創(chuàng)建顏色圖
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)
#添加格網(wǎng)
g3=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g3.xformatter = LONGITUDE_FORMATTER  # x軸設(shè)為經(jīng)度的格式
g3.yformatter = LATITUDE_FORMATTER  # y軸設(shè)為緯度的格式
#設(shè)置經(jīng)緯度網(wǎng)格的間隔
g3.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g3.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改經(jīng)緯度字體大小 
g3.xlabel_style={'size':10,'weight': 'bold'}                            
g3.ylabel_style={'size':10,'weight': 'bold'}
#添加標(biāo)題
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()

6.添加圖例

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)
#添加格網(wǎng)
g3=ax.gridlines(linestyle='--',draw_labels=True,linewidth=1,alpha=0.5)
g3.xformatter = LONGITUDE_FORMATTER  # x軸設(shè)為經(jīng)度的格式
g3.yformatter = LATITUDE_FORMATTER  # y軸設(shè)為緯度的格式
#設(shè)置經(jīng)緯度網(wǎng)格的間隔
g3.xlocator = mticker.FixedLocator(np.arange(-180, 180, 60))
g3.ylocator = mticker.FixedLocator(np.arange(-90, 90, 30))
#修改經(jīng)緯度字體大小 
g3.xlabel_style={'size':10,'weight': 'bold'}                            
g3.ylabel_style={'size':10,'weight': 'bold'}
#添加標(biāo)題
ax.set_title('2018 CO$_2$ Emissions',
              fontsize=15,
              fontweight ='bold',
              pad=20)
#設(shè)置圖例
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()

7.參考資料

Crippa, M., Guizzardi, D., Schaaf, E., Solazzo, E., Muntean, M., Monforti-Ferrario, F., Olivier, J.G.J., Vignati, E.: Fossil CO2 and GHG emissions of all world countries — 2021 Report, in prep.

Crippa, M., Solazzo, E., Huang, G., Guizzardi, D., Koffi, E., Muntean, M., Schieberle, C., Friedrich, R. and Janssens-Maenhout, G.: High resolution temporal profiles in the Emissions Database for Global Atmospheric Research. Sci Data 7, 121 (2020). doi:10.1038/s41597–020–0462–2.

IEA (2019) World Energy Balances, www.iea.org/data-and-statistics, All rights reserved, as modified by Joint Research Centre, European Commission.

Jalkanen, J. P., Johansson, L., Kukkonen, J., Brink, A., Kalli, J., & Stipa, T. (2012). Extension of an assessment model of ship traffic exhaust emissions for particulate matter and carbon monoxide. Atmospheric Chemistry and Physics, 12(5), 2641–2659. doi:10.5194/acp-12–2641–2012

Johansson, L., Jalkanen, J.-P., & Kukkonen, J. (2017). Global assessment of shipping emissions in 2015 on a high spatial and temporal resolution. Atmospheric Environment, 167, 403–415. doi:10.1016/j.atmosenv.2017.08.042

https:///visualising-the-worlds-carbon-dioxide-emissions-with-python-e9149492e820

費(fèi)老師geopandas教程、云臺(tái)cartopy教程

https://github.com/hhuangwx/cmaps

Cartopy. Met Office. git@github.com:SciTools/cartopy.git. 2015-02-18. 7b2242e.

最后，您可以在公眾號(hào)輸入“二氧化碳”，可以得到測(cè)試數(shù)據(jù)。

希望關(guān)注的朋友們轉(zhuǎn)發(fā)，也可以將微信公眾號(hào)設(shè)置為星標(biāo)，下次我們更新發(fā)送的推文您就可以及時(shí)的查看到了。如果不方便轉(zhuǎn)發(fā)或者設(shè)置星標(biāo)，那就請(qǐng)點(diǎn)個(gè)贊或者再看。

前端百思說(shuō)

不喜歡

AI小白筆記

不喜歡

AI小白筆記

不喜歡