编程分享|利用Python可视化全球二氧化碳排放
原创 GuaGua 农业与土地遥感 2022-05-19 08:30 发表于美国
收录于合集
#编程分享9
#遥感美图1
#Python1
背景
近几年碳的问题一直很火，那么关于碳数据的收集就比较重要，刚好前不久见到一个全球实时碳数据的网站，里面包含碳的数据，网站包含一些主要国家碳排放数据，碳排放数据涉及到的来源有： “电力”、“地面运输”、“工业”、“居民消费”、“国内航空”、“国际航空”。
网站链接：https:///
此网站不仅数据好，而且图也很漂亮，刚打开首页，就被图吸引住了。
那么我们如何绘制地图呢？刚好有个国外大神Adam Symington画了很多好看的地图，其中就有二氧化碳排放的地图。我们借鉴模仿大神的思路，选择自己想要表达的效果。
Adam Symington作品
首先是获取数据，数据来源同大神作品的数据，来源于另一个碳数据网站：https://edgar.jrc.ec./ 此网站有许多数据集可供您选择，不仅仅有碳排放数据。
在本文中，我们只研究 2018 年全球的二氧化碳排放量，其他数据不进行研究，感兴趣的可以去原网站探索。数据的描述和原理见原网站，本文不再做赘述，直接开始对数据进行可视化。
1.读取二氧化碳排放数据
#!/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
2.可视化画图
#首先进行简单可视化
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()
3.对图形进行优化一下
#对数颜色映射的一个小例子，帮助大家理解下面颜色映射，为后续做好准备
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()
4.对图形进行投影
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
| lat | lon | emission | geometry |
|---|
| 0 | 88.1 | -50.7 | 11.0301 | POINT (-50.70000 88.10000) |
|---|
| 1 | 88.1 | -50.6 | 15.7066 | POINT (-50.60000 88.10000) |
|---|
| 2 | 88.1 | -50.5 | 19.1943 | POINT (-50.50000 88.10000) |
|---|
| 3 | 88.1 | -50.4 | 22.0621 | POINT (-50.40000 88.10000) |
|---|
| 4 | 88.1 | -50.3 | 24.5536 | POINT (-50.30000 88.10000) |
|---|
| ... | ... | ... | ... | ... |
|---|
| 2191880 | -77.8 | 166.5 | 11459.2000 | POINT (166.50000 -77.80000) |
|---|
| 2191881 | -77.8 | 166.6 | 694.4980 | POINT (166.60000 -77.80000) |
|---|
| 2191882 | -77.9 | 166.3 | 2133.1000 | POINT (166.30000 -77.90000) |
|---|
| 2191883 | -77.9 | 166.4 | 1389.0000 | POINT (166.40000 -77.90000) |
|---|
| 2191884 | -77.9 | 166.5 | 7887.5100 | POINT (166.50000 -77.90000) |
|---|
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首先得定义坐标系
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()
#使用Robinson投影，以中国为中心
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)
#添加格网
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.show()
5.色彩映射
#发现前面颜色对比不够明显，我们重新设置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()
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)
#添加格网
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()
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
费老师geopandas教程、云台cartopy教程
https://github.com/hhuangwx/cmaps
Cartopy. Met Office. git@github.com:SciTools/cartopy.git. 2015-02-18. 7b2242e.
最后，您可以在公众号输入“二氧化碳”，可以得到测试数据。
希望关注的朋友们转发，也可以将微信公众号设置为星标，下次我们更新发送的推文您就可以及时的查看到了。如果不方便转发或者设置星标，那就请点个赞或者再看。
农业与土地遥感
分享农业与土地遥感领域，包括可视化、编程、农业遥感及其交叉学科领域的最新科研进展。本公众号由中国农业大学土地科学与技术学院农业与土地遥感团队维护与更新。
18篇原创内容
公众号
