Python遥感可视化 — Basemap将地面观测站点进行空间插值可视化

欢迎关注博主的微信公众号：“智能遥感”。

该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。

本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

“本节主要内容是将地面观测站点数据在地图上显示，并且对站点数据进行空间插值，从而可以将点数据转化为栅格数据，栅格数据具有较好的空间连续性。”

今天的遥感之美封面图—白雪下的麦肯齐河系统。乍一看，有一种山回路转不见君，雪上空留马行处的景象。麦肯齐河系统是加拿大最大的流域，也是世界上第十大流域。这条河从加拿大落基山脉的哥伦比亚冰原到达北冰洋，全长4,200公里。

2016年11月7日由Landsat 8上OLI传感器拍摄的图像，该图显示了Mackenzie River Delta的一部分 - 包括沿着河道东部的这条冰冻高速公路。由白色冰雪覆盖的水道在绿色松树覆盖的土地上脱颖而出。

PM2.5可视化

在平时遥感数据可视化时，等经纬度投影是我们最长使用的一种方法，也是最简单的投影。等经纬度投影又称等距圆柱投影，是假想球面与圆筒面相切于赤道，赤道为没有变形的线。经纬线网格，同一般正轴圆柱投影，经纬线投影成两组相互垂直的平行直线。其特性是：保持经距和纬距相等，经纬线成正方形网格；沿经线方向无长度变形；角度和面积等变形线与纬线平行，变形值由赤道向高纬逐渐增大。该投影适合于低纬地区制图。

首先来看一下Basemap官网上的参考代码，这是最简单的显示全球地图的方式（图1）：大陆廓线（未添加国家行政边界）以及陆地（橘黄色）、海洋（浅绿色）、内陆湖泊（蓝色）的颜色填充，让人更加直观地了解各大洲的分布情况。

代码实现：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2019/1/11 14:09'from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as pltmaps = Basemap(projection='cyl', lat_0=0, lon_0=0)
# 首先给地球涂上蓝色
maps.drawmapboundary(fill_color='aqua')
# 再给大陆涂上橘黄色,给江河湖泊涂上大海一样的颜色
maps.fillcontinents(color='coral', lake_color='blue')
maps.drawcoastlines()
plt.show()
# plt.savefig('test.png')

结果图：

图1 等经纬度投影下的世界地图

在上述基础之上，我们将地面观测站点数据叠加到底图上，来分析数据的时空变化。下面以2018年01月01日中国地区PM2.5观测数据为例，图2展示全国PM2.5站点空间分布，我们发现中国东部地区观测站点数量加多，并且点相对较大，点的大小对应PM2.5数值的高低，说明东部地区这一天的大气污染较西部地区较为严重。

图2 全国PM2.5站点分布，点越大说明PM2.5浓度越高

但是上图不能够很好地反应PM2.5的空间变化，从视觉效果来看，颜色较为单一。因此，我们给每个站点添加颜色，颜色越红，表示PM2.5浓度越大，也就是空气质量越差（见下图）。点的颜色和大小可以同时表示数值大小，更加增强了数据的可视化效果。

图3 全国PM2.5站点分布，颜色越红且点越大说明PM2.5浓度越高

代码实现：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2019/1/11 14:49'from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np# 用来正常显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
# 用来正常显示负号
plt.rcParams['axes.unicode_minus'] = False
# 获取PM2.5数据
df = pd.read_excel(r'D:\data\20180101PM25-CHINA.xlsx')
# 剔除无效值NAN
df = df.dropna(axis=0, how='any')lat = np.array(df["lat"][:])  # 获取维度之维度值
lon = np.array(df["lon"][:])  # 获取经度值
PM25 = np.array(df["PM25"][:], dtype=float)
# 画图
fig = plt.figure(figsize=(16, 9))
plt.rc('font', size=15, weight='bold')
ax = fig.add_subplot(111)
# 添加标题,PM2.5下标设置
plt.title(u'2018年01月01日中国地区$\mathrm{PM}_{2.5}$质量浓度分布', size=25, weight='bold')
# 创建底图,等经纬度投影
mp = Basemap(llcrnrlon=73., llcrnrlat=17.,urcrnrlon=135., urcrnrlat=55,projection='cyl', resolution='h')
# 添加海岸线
mp.drawcoastlines()
# 添加国家行政边界
mp.drawcountries()
# 设置colorbar中颜色间隔个数
levels = np.linspace(0, np.max(PM25), 20)
# cf = mp.scatter(lon, lat, PM25, marker='o', color='r')
# 设置颜色表示数值大小
cf = mp.scatter(lon, lat, PM25, c=PM25, cmap='jet', alpha=0.75)
# 设置上下标以及单位的希腊字母
cbar = mp.colorbar(cf, location='right', format='%d', size=0.3,ticks=np.linspace(0, np.max(PM25), 10),label='$\mathrm{PM}_{2.5}$($\mu$g/$\mathrm{m}^{3}$)')
plt.show()

地面站点数据观测精度相对较高，但是不能够覆盖所有的空间范围。为了研究大范围连续空间的特征变化，通常需要对站点数据进行空间插值，例如克里金插值、样条插值、反距离加权。下面以三种插值为例（最邻近插值、线性插值和三次样条插值）。如果喜欢克里金插值的小伙伴，可以pykrige模块去实现，本节不作讲解。

图4 最邻近插值后，全国PM2.5空间分布，颜色越红说明PM2.5浓度越高

图5 线性插值后，全国PM2.5空间分布，颜色越红说明PM2.5浓度越高

图6 三次样条插值后，全国PM2.5空间分布，颜色越红说明PM2.5浓度越高

说明：三次样条插值后的效果其实分布也是连续的，与图5很相似，但是在Basemap上显示就出现破碎现象，目前还不知道什么原因。

代码实现：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2019/1/11 14:49'from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from scipy.interpolate import griddata
import cmaps# 用来正常显示中文
plt.rcParams['font.sans-serif'] = ['SimHei']
# 用来正常显示负号
plt.rcParams['axes.unicode_minus'] = False
# 获取PM2.5数据
df = pd.read_excel(r'D:\data\20180101PM25-CHINA.xlsx')
# 剔除无效值NAN
df = df.dropna(axis=0, how='any')
# 获取纬度
lat = np.array(df["lat"][:])
# 获取经度
lon = np.array(df["lon"][:])
# 获取PM2.5
PM25 = np.array(df["PM25"][:])
# 创建格网
grid_x, grid_y = np.mgrid[73.56:135.04, 18.2:53.56]
# 插值方法：'nearest', 'linear', 'cubic'
grid_z = griddata((lon, lat), PM25, (grid_x, grid_y), method='cubic')# 画图
fig = plt.figure(figsize=(16, 9))
plt.rc('font', size=15, weight='bold')
ax = fig.add_subplot(111)
# 添加标题,PM2.5下标设置
plt.title(u'2018年01月01日中国地区$\mathrm{PM}_{2.5}$质量浓度分布', size=25, weight='bold')
# 创建底图,等经纬度投影
mp = Basemap(llcrnrlon=73., llcrnrlat=17.,urcrnrlon=135., urcrnrlat=55,projection='cyl', resolution='h')
# 添加海岸线
mp.drawcoastlines()
# 添加国家行政边界
mp.drawcountries()
# 设置colorbar中颜色间隔个数
levels = np.linspace(0, np.max(PM25), 20)
cf = mp.contourf(grid_x, grid_y, grid_z, levels=levels, cmap=cmaps.GMT_panoply)
cbar = mp.colorbar(cf, location='right', format='%d', size=0.3,ticks=np.linspace(0, np.max(PM25), 10),label='$\mathrm{PM}_{2.5}$($\mu$g/$\mathrm{m}^{3}$)')
plt.show()

Python遥感可视化 — Basemap将地面观测站点进行空间插值可视化相关推荐

Python遥感可视化 — Basemap对遥感数据可视化
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
Python遥感可视化 — folium模块展示热力图
欢迎关注博主的微信公众号:"智能遥感". 该公众号将为您奉上Python地学分析.爬虫.数据分析.Web开发.机器学习.深度学习等热门源代码. 本人的GitHub代码资料主页(持续 ...
python 克里金空间插值_Python克里金(Kriging)插值计算及可视化绘制
前面两篇推文我们分别介绍了使用Python和R进行IDW(反距离加权法) 插值的计算及结果的可视化过程,详细内容可见如下: 本期推文,我们将介绍如何使用Python进行克里金(Kriging)插值计算 ...
使用python爬取BOSS直聘岗位数据并做可视化（Boss直聘对网页做了一些修改，现在的代码已经不能用了）
使用python爬取BOSS直聘岗位数据并做可视化结果展示首页岗位信息岗位详情薪资表学历需求公司排名岗位关键词福利关键词代码展示爬虫代码一.导入库二.爬取数据 1.爬取数据代 ...
Python使用matplotlib可视化多个时间序列数据、在同一个可视化图像中可视化多个时间序列数据（Multiple Time Series）
Python使用matplotlib可视化多个时间序列数据.在同一个可视化图像中可视化多个时间序列数据(Multiple Time Series) 目录
python使用matplotlib绘制一条正弦曲线（plot函数可视化sine plot）
python使用matplotlib绘制一条正弦曲线(plot函数可视化sine plot) 目录 python使用matplotlib绘制一条正弦曲线(plot函数可视化sine plot) #导入 ...
Python在Seaborn中手动指定调色板颜色进行数据可视化颜色自定义实战（Manually Specify Palette Colors in Seaborn）
Python在Seaborn中手动指定调色板颜色进行数据可视化颜色自定义实战(Manually Specify Palette Colors in Seaborn) 目录
python使用matplotlib可视化线图（line plot）、为可视化图像添加双Y轴、分别可视化不同范围的数据（double y axis in matplotlib）
python使用matplotlib可视化线图(line plot).为可视化图像添加双Y轴.分别可视化不同范围的数据(double y axis in matplotlib) 目录
seaborn可视化条形图并按照升序排序条形图进行可视化：Sort Bars in Barplot in Ascending Order in Python
seaborn可视化条形图并按照升序排序条形图进行可视化:Sort Bars in Barplot in Ascending Order in Python 目录

Python遥感可视化 — Basemap将地面观测站点进行空间插值可视化

Python遥感可视化 — Basemap将地面观测站点进行空间插值可视化相关推荐

最新文章

热门文章