引言
好久好久没有更新了，来杭州的这几天真的是快把人烤熟了，好在学校的空调给力，连厕所都是中央空调，有钱真的会玩呀......
继上一篇博客《Python处理GPM(IMERG/GSMaP)卫星降水数据》 下面说说如何利用python读取GPM IMERG降水数据并绘制降水分布图
- Basemap绘图工具包,自行下载安装：python工具包下载地址：
https://www.lfd.uci.edu/~gohlke/pythonlibs/#pymssql - 作者以IMERG 日尺度降水数据为例子，该示例数据可从NASA网站上下载：
https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/2000/10/3B-DAY.MS.MRG.3IMERG.20001026-S000000-E235959.V06.nc4 -
代码示例
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @File : gpm_mapping.py
# @Author : zengsk in HHU
# @Time : 2019/7/21 13:40
import h5py
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from netCDF4 import Dataset as ncdataset
from mpl_toolkits.basemap import Basemap
####### 绘制GPM日降水分布图 ###########
'''
自定义色带(可替换 map.contourf()函数的cmap参数)
ps：我自己定义的实在难看，还是用自带的吧^-^
'''
def colormap():
ColorList = ['#DDDDFF','#7D7DFF','#0000C6','#000079','#CEFFCE',
'#28FF28','#007500','#FFFF93','#8C8C00','#FFB5B5',
'#FF0000','#CE0000','#750000']
return mpl.colors.LinearSegmentedColormap.from_list('cmap', ColorList, 256)
filename = r'../assets/3B-DAY.MS.MRG.3IMERG.20001026-S000000-E235959.V06.nc4'
# 这部分是读取小时分辨率数据的代码
# fHander = h5py.File(filename, mode='r')
# precip = fHander['/Grid/precipitationCal'][:]
# precip = np.transpose(precip * 0.5)
# lat = fHander['/Grid/lat'][:]
# lon = fHander['/Grid/lon'][:]
dst = ncdataset(filename)
precip = dst.variables['precipitationCal'][:]
precip = precip.reshape(3600, 1800)
precip = np.transpose(precip)
precip[np.isnan(precip)] = -999
lat = dst.variables['lat'][:]
lon = dst.variables['lon'][:]
# Plot the figure, define the geographic bounds
fig = plt.figure(3, dpi=300)
map = Basemap()
# Draw coastlines, state and country boundaries, edge of map.
map.drawcoastlines(linewidth=1.2, linestyle='solid', color='k', antialiased=3)
map.drawcountries(linewidth=0.8, linestyle='solid', color='k', antialiased=3)
# map.bluemarble(scale=3) # 加载蓝色大理石背景
# map.shadedrelief(scale=2) # 加载地形阴影
# map.etopo(scale=3) # 加载地形包括海洋
# Draw filled contours.
clevs = np.arange(0, 100, 10)
# color scale changing
# clevs = [0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
# Define the latitude and longitude data
x, y = np.float32(np.meshgrid(lon, lat))
# Mask the values less than 0 because there is no data to plot.
masked_array = np.ma.masked_where(precip < 0, precip)
# Plot every masked value as white
# colormap
cmap = plt.cm.get_cmap('jet')
cmap.set_bad('w', .8)
# Plot the data
cs = map.contourf(x, y, precip, clevs, cmap=cmap, latlon=True)
parallels = np.arange(-60., 61, 20.)
map.drawparallels(parallels, labels=[True, False, True, False])
meridians = np.arange(-180., 180., 60.)
map.drawmeridians(meridians, labels=[False, False, False, True])
# 设置标题和字体
plt.title('GPM IMERG Daily Rainfall( mm )')
font = {'weight': 'bold', 'size': 6}
plt.rc('font', **font)
# 加载色带
cbar = map.colorbar(cs, location='right', pad="5%")
cbar.set_label('mm/daily')
plt.savefig('../output/gpm.jpg', dpi=300) # change to your directory
plt.show()
print(" ###### 数据处理完成 ###### ")
结果展示：
image.png
