Python | GPM CMB数据求降水月平均

绘图目标如标题

数据格式为hdf5

1. 数据处理

import h5py
from importlib.resources import path
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import time

time_start = time.time()

#数据处理
def data_clean(path,extent):
    with h5py.File(str(path), 'r') as f:
        lon = pd.DataFrame(f['NS/Longitude'][:])
        lat = pd.DataFrame(f['NS/Latitude'][:])
        surfacePrecipitation = pd.DataFrame(f['NS/surfPrecipTotRate'][:])
    surfacePrecipitation.replace(-9999.9,np.nan)
    time1 = str(path)[38:64]
    lonmin, lonmax, latmin, latmax = extent
    mask = (lon >= lonmin) & (lon <= lonmax) & (lat >= latmin) & (lat <= latmax)
    x = pd.DataFrame(mask)
    lon_s = lon[x].dropna(axis=0,how='all').dropna(axis=1,how='all') #lon[x]后数据size不变，dropna后裁剪到需要的范围，size变小
    lat_s = lat[x].dropna(axis=0,how='all').dropna(axis=1,how='all')
    sp_s = surfacePrecipitation[x].dropna(axis=0,how='all').dropna(axis=1,how='all')
    return lon_s,lat_s,sp_s,time1

def sum1(a,b):
    s=0
    if (pd.isnull(a) == True):
        s=b 
    elif(pd.isnull(b) == True):
        s=a
    else:
        s=a+b
    return s

def sum2(a,b):
    p1 ,q1= a.shape
    p2,q2 = b.shape
    if ((p1 !=p2) | (q1 != q2)):
        print('Wrong!')
    else:
        s=np.zeros([p1,q1])
        for i in range(0,p1):
            for j in range(0,q1):
                if (pd.isnull(a[i][j]) == True):
                    s[i][j]=b[i][j] 
                elif(pd.isnull(b[i][j]) == True):
                    s[i][j]=a[i][j] 
                else:
                    s[i][j]=b[i][j] +a[i][j]

    return s

#区域设置
extent =[106,111,26,31]
lonn1,lonx1,latn1,latx1 = extent
dx = lonx1 - lonn1
dy = latx1 - latn1
#分辨率设置
dx_small = 0.25 
dy_small = 0.25
n1 = int(dx/dx_small) #水平小格数
n2 = int(dy/dy_small) #竖直小格数

path = r"D:\桌面\CMB样本"
files= os.listdir(path)
time_total = []  #记录每个文件代表的扫描时间

sp_all = np.full([n2,n1], np.nan)  #二维数组，记录所有文件每个小格里的总降水
scan_times = np.zeros([n2,n1]) #二维数组，记录所有文件每个小格里的footprints
for filename in files:
    sp_one = np.full([n2,n1], np.nan) #二维数组，记录单个文件每个小格里的总降水
    number = np.zeros([n2,n1]) #二维数组，记录单个文件每个小格里的footprints
    lon,lat,sp,time1 = data_clean(os.path.join(path,filename),extent) #获得裁剪后符合区域的数据 ，可以使循环变短
    imax,jmax = lon.shape
    if (imax!=0): #裁剪后有数据，数据中含nan（没扫描到的地方）
        for i in range(0,imax):
            for j in range(0,jmax):
                if(pd.isnull(lat.iloc[i:(i+1),j:(j+1)].values) != True): #lat.iloc[i:(i+1),j:(j+1)]是DataFrame格式数据的索引方式，得到的是一个纬度
                    #pd.isnull()！=True说明这个数不是nan，是扫描经过的点，下面判断这个点在哪个格子里
                    a = lat.iloc[i:(i+1),j:(j+1)].values
                    b = lon.iloc[i:(i+1),j:(j+1)].values
                    e = n2-1-int(((a-latn1)//dy_small))
                    f = int(((b-lonn1)//dx_small))
                    sp_one[e][f] = sum1(float(sp.iloc[i:(i+1),j:(j+1)].values),sp_one[e][f]) #sum1是自己写的函数，因为相加时会有nan参与的情况，但dataframe的计算不支持
                    number[e][f] = number[e][f] +1     #小格里的footprints数量+1        
    time_total.append(time1) #记录扫描时间
    for i in range(0,n2):
        for j in range(0,n1):
            if number[i][j] !=0:
                sp_one[i][j] = sp_one[i][j]/number[i][j]
    sp_all = sum2(sp_all,sp_one) #sum2也是自己写的函数，因为相加时会有nan参与的情况，但dataframe的计算不支持
    
    scan_times[~pd.isnull(sp_one)] +=1 #spv_one不是nan的数的位置上scan_times+1，这次扫描这个格子里有数据（不是nan）

np.savetxt(r"D:\桌面\20140703.txt", sp_all, delimiter ="\t")
np.savetxt(r"D:\桌面\201407scantimes.txt", scan_times, delimiter ="\t")

time_end = time.time()  # 记录结束时间
time_sum = time_end - time_start  # 计算的时间差为程序的执行时间，单位为秒/s
print('运行时间:',time_sum)

2. 绘图

import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter


def make_Rr_cmap(levels):
    '''制作降水的colormap.'''
    nbin = len(levels) - 1
    cmap = cm.get_cmap('jet', nbin)
    norm = mcolors.BoundaryNorm(levels, nbin)
    return cmap, norm



extents =[106,111,26,31]
lonn1,lonx1,latn1,latx1 = extents
arr = np.loadtxt(r"D:\20140703.txt")


lon1 = np.arange(lonn1,lonx1,0.25)
lat1 = np.arange(latn1,latx1,0.25)
lon,lat = np.meshgrid(lon1,lat1)

levels_Rr = [0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 15.0, 20.0, 25]
#levels_Rr = [0.0,0.25, 0.5,0.75, 1.0, 1.25,1.5,1.75,2.0]
cmap_Rr, norm_Rr = make_Rr_cmap(levels_Rr)

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=proj)
ax.set(xlim=(extents[0], extents[1]),ylim=(extents[-2], extents[-1]))
ax.tick_params(labelsize='large')
im = plt.imshow(arr,cmap=cmap_Rr, norm=norm_Rr,extent=extents,zorder=3,alpha=0.8)
china = Reader(r"D:\bou2_4l.dbf").geometries()
ax.add_geometries(china, ccrs.PlateCarree(),
                  facecolor='none', edgecolor='black', zorder=2,linewidth=1.5) 
river = Reader(r"D:\1级河流.shp").geometries()
ax.add_geometries(river, ccrs.PlateCarree(), facecolor='none', edgecolor='RoyalBlue', linewidth=1)

cb = fig.colorbar(im, ax=ax, extend='both')
cb.set_label('Rain Rate (mm/hr)', fontsize=15)
cb.ax.tick_params(labelsize=15)
ax.set_xticks(np.arange(extents[0], extents[1]+1, 1), crs=proj)
ax.set_yticks(np.arange(extents[-2], extents[-1]+1, 1), crs=proj)

ax.xaxis.set_major_formatter(LongitudeFormatter())#zero_direction_label=False))
ax.yaxis.set_major_formatter(LatitudeFormatter())
miloc = plt.MultipleLocator(0.25) #0.5意思是以0.5°为间距画线
ax.xaxis.set_minor_locator(miloc)
ax.grid(axis='x', which='both') #both可替换为major或minor，画的线不同，可自行尝试
 
miloc = plt.MultipleLocator(0.25)
ax.yaxis.set_minor_locator(miloc)
ax.grid(axis='y', which='both')

ax.set_title('20140703 SurfPrecipTotRate', fontsize=20,pad=20) #on Average
plt.tick_params(labelsize=15)
plt.show()

出图：