Python遥感可视化 — Basemap对遥感数据可视化

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

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

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

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

---

“小伙伴们，好久不见！本节接着上一节的内容继续分享遥感可视化功能，主要通过对欧洲再分析资料栅格数据进行可视化。”

今天的遥感之美 — 罗平县，中国的“黄金土地”。罗平县隶属于云南省曲靖市，地处东经103°57′～104°43′、北纬24°31′～25°25′之间，除南部八大河一带属南亚热带气候外，其余为高原季风气候。每到春天来临之时，正是油菜花绽放的最佳时期。

每年的2月中旬，卫星图像中的黄色区域在我国南方罗平县的山谷中蔓延，这便是那可爱的油菜花短暂的绽放。

使用中分辨率成像光谱仪（MODIS）的粗分辨率（每个像元为250米），效果是油菜花绽放的山谷比它们周围的森林山丘和山脊的深绿色稍微淡了些。

凭借Landsat 8上的Operational Land Imager（OLI）更精细的分辨率（每像元为30米），效果更加清晰。广阔的黄色田野和牙齿形状的喀斯特山丘的绿色之间的界限是显而易见的。

Landsat 8于2017年2月19日获得了上图，通常这是油菜花处于一年中的高峰期。到了三月中旬，花朵逐渐消失，绿色重新凸起。

---

在上一次的Basemap可视化分享中，主要将地面站点以及空间插值后的数据进行可视化。这一节，主要对栅格数据直接可视化。在日常遥感数据可视化中，遥感栅格、在分析资料等都是常见的数据，下面以欧洲中期天气预报中心（ECMWF）的再分析资料为例，来实现Basemap栅格可视化功能。

ECMWF下载地址：https://apps.ecmwf.int/datasets/，关于如何运用Python对再分析资料批量下载，网上博客已有相关介绍，这个在后面数据爬虫部分，我也会详细介绍，就请大家耐心等待，后面还有很多干货！

ECMWF再分析资料包含温、压、湿、风、流场、数值预报产品等,一般主要看他的数值预报,如流场、温度、气压场、降水等。除了这些常见气象数据，还有云、气溶胶、痕量气体、辐射、地表径流等。

代码实现：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2019/3/10 21:22'

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import cmaps

meteo_file = r'D:\data\ECMWF.nc'
ds = Dataset(meteo_file, mode='r')

# 获取每个变量的值
lons = ds.variables['longitude'][:]
lats = ds.variables['latitude'][:]

# surface_air_pressure
sp = ds.variables['sp'][:]
sp_units = ds.variables['sp'].units
scale_factor = ds.variables['sp'].scale_factor
add_offset = ds.variables['sp'].add_offset
sp = scale_factor * sp + add_offset
# 2 metre temperature
t2m = ds.variables['t2m'][:]
t2m_units = ds.variables['t2m'].units
scale_factor = ds.variables['t2m'].scale_factor
add_offset = ds.variables['t2m'].add_offset
t2m = scale_factor * t2m + add_offset
# Total column ozone
tco3 = ds.variables['tco3'][:]
tco3_units = ds.variables['tco3'].units
scale_factor = ds.variables['tco3'].scale_factor
add_offset = ds.variables['tco3'].add_offset
tco3 = scale_factor * tco3 + add_offset

# 经纬度平均值
lon_0 = lons.mean()
lat_0 = lats.mean()

# 画图大小设置
fig = plt.figure(figsize=(16, 9))
plt.rc('font', size=15, weight='bold')
ax = fig.add_subplot(111)

m = Basemap(lat_0=lat_0, lon_0=lon_0)
lon, lat = np.meshgrid(lons, lats)
xi, yi = m(lon, lat)


# 这里数据时间是UTC 00:00，2018年1月的日平均数据，只展示1月1号的数据
sp_01 = sp[0:1, :, :]
t2m_01 = t2m[0:1, :, :]
tco3_01 = tco3[0:1, :, :]

levels = m.pcolor(xi, yi, np.squeeze(tco3_01), cmap=cmaps.GMT_panoply)

# 添加格网与绘制经纬线
m.drawparallels(np.arange(-90., 91., 20.), labels=[1, 0, 0, 0], fontsize=15)
m.drawmeridians(np.arange(-180., 181., 40.), labels=[0, 0, 0, 1], fontsize=15)

# 添加海岸线,省州边界以及国家行政边界
m.drawcoastlines()
m.drawstates()
m.drawcountries()

# 添加colorbar
cbar = m.colorbar(levels, location='bottom', pad="10%")
cbar.set_label(tco3_units, fontsize=15, weight='bold')

# 添加图的标题
plt.title('Total column ozone')
plt.show()

ds.close()

结果图：

2018年01月01日近地面大气压强（单位：Pa，分辨率：1°×1°）

2018年01月01日近2米处温度（单位：K，分辨率：1°×1°）

2018年01月01日臭氧总柱状浓度（单位：kg/m-2，分辨率：1°×1°）