[Python知识库]（本篇仅用作本人记录写过的代码）GEE（python）批量降水产品转数组，并下载，以及分析降水数据频次强度

写这篇博客的心情是郁闷的，三个月以来因为研究思路和经验不足的问题，一直在推翻重做。这一次也毫不例外…跑了半个月的代码，又付诸东流了。
本篇代码基于 较新的定义方法 小时极端降水的事件的定义：一小时及以上降水超过80%百分位值，连续8小时不超过80%百分值视为该事件结束，并且当日降水量达到90%百分位值。
可惜最终出图效果差，被pass了。接下来记录，分析数据所写的代码：

批量降水产品转数组

#Edited by Xinglu Cheng 2022.01.09
#将每年的半小时降水数据转为三维数组，分别保存到本地
import ee
import os
import threading
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
from ipygee import Map
from IPython.display import Image
import sys
sys.setrecursionlimit(10000)#修改递归深度

import geemap
import geemap.colormaps as cm
gee_Map = geemap.Map(center=[36.56, 116.98], zoom=6)
gee_Map.add_basemap('Gaode.Normal')

#加入研究区------------------------------------------------------------------------------------------------
HHH = ee.FeatureCollection("users/2210902126/GPM/HHH_Area")

#逐天获取图像，转数组    
def toDailyArray(year,DailyArray,band):
    for i in range(1,13):
        print(i)
        if i == 1 or i == 3 or i == 5 or i == 7 or i == 8 or i == 10 or i == 12: day = 31
        if i == 4 or i == 6 or i == 9 or i == 11: day = 30
        if i == 2 and year % 4 != 0 or i == 2 and year % 4 == 0 and year % 400 != 0: day = 28
        if i == 2 and year % 4 == 0 and year % 100 != 0 or i == 2 and year % 4 == 0 and year % 400 == 0: day = 29
         # 判断一个月的天数

        for d in range(1, day+1):
            for t in range(0,24):
                if t<10:
                    if i >= 10 and day >= 10: extend = ee.Date(str(year) + '-' + str(i) + '-' + str(d)+'T0'+str(t)+':00:00')
                    if i >= 10 and day < 10: extend = ee.Date(str(year) + '-' + str(i) + '-0' + str(d)+'T0'+str(t)+':00:00')
                    if i < 10 and day >= 10: extend = ee.Date(str(year) + '-0' + str(i) + '-' + str(d)+'T0'+str(t)+':00:00')
                    if i < 10 and day < 10: extend = ee.Date(str(year) + '-0' + str(i) + '-0' + str(d)+'T0'+str(t)+':00:00')
                if t>=10:
                    if i >= 10 and day >= 10: extend = ee.Date(str(year) + '-' + str(i) + '-' + str(d)+'T'+str(t)+':00:00')
                    if i >= 10 and day < 10: extend = ee.Date(str(year) + '-' + str(i) + '-0' + str(d)+'T'+str(t)+':00:00')
                    if i < 10 and day >= 10: extend = ee.Date(str(year) + '-0' + str(i) + '-' + str(d)+'T'+str(t)+':00:00')
                    if i < 10 and day < 10: extend = ee.Date(str(year) + '-0' + str(i) + '-0' + str(d)+'T'+str(t)+':00:00')
            # 判断月份是否是两位数，两位数不需要加0,同理并判断小时
                if d==8 and t==3:
                    band=band
                else:
                    image=ee.ImageCollection("NASA/GPM_L3/IMERG_V06").filter(ee.Filter.date(extend.getRange('hour'))).select('precipitationCal').sum().multiply(0.5).clip(HHH).reproject(crs='EPSG:4326',scale=11132)
                    #获取图像，裁切重采样

                    if band==0:
                        numArray=geemap.ee_to_numpy(image,region=HHH,default_value=-999)
                        DailyArray=numArray.copy()
                    else:
                        numArray=geemap.ee_to_numpy(image,region=HHH,default_value=-999).squeeze()#图像转数组
                        DailyArray=np.insert(DailyArray,band,numArray,axis=2)#将数组存入
                    band+=1
    return DailyArray,band

    
#逐年获取数据
newpath=os.path.expanduser('~/gee_data')
year_start=2001
year_end=2022
bands=0
for i in range(year_start,year_end):
    if i == year_start:
        DailyArray,bands=toDailyArray(i,bands,bands)
    else:
        DailyArray,bands=toDailyArray(i,DailyArray,bands)
    print(i)
    # 存储三维数组
    newpath = os.path.join(newpath, "Hour_data"+str(i)+"_12.npy")
    np.save(file=newpath, arr=DailyArray)
print(bands)

分析降水数据频次强度

#Edited by Xinglu Cheng 2022.01.27
#读取小时数据，并统计极端降水特征
import ee
import os
import threading
import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal
import sys
sys.setrecursionlimit(10000)#修改递归深度

#通过小时降水数据，获取日降水百分值90%--------------------------------------------------------------------
def GetDailyPre(array,band):
    for t in range(23,band,24):
        if t==23:
            daily_array=np.array([np.sum(array[t-23:t])])
        else:
            daily_array=np.append(daily_array,np.sum(array[t-23:t]))
    width = round(daily_array.shape[0]*0.9)#四舍五入找到第90%的位置
    Threshold = np.sort(daily_array)[width]
    return Threshold

#判定极端降水事件,并统计时长-----------------------------------------------------------------------------
def DefineEx(array,band):
    num = 0
    numm = 0
    one_pre = 0 #单次降水量
    pre_array = np.zero((1))#极端降水量数组
    dura_data = 0 #总极端降水量
    fre_data = 0 #频次
    num_data = 0 #总时长
    dailyThre=GetDailyPre(array,band)#获取日降水阈值    
    for n in range(band):
        width = round(array.shape[0]*0.8)#四舍五入找到第80%的位置
        Threshold = np.sort(array)[width]
        #判断单次极端降水事件
        if array[n] >= Threshold:
            dura_data += array[n]
            one_pre  += array[n]
            num += 1
            numm = 0
        if num>0 and array[n] < Threshold and numm<7:
            numm += 1
        if num>0 and numm == 7 and array[n] < Threshold:
            if np.sum(array[(n-8-round(num/2))-12:(n-8-round(num/2))+11])>=dailyThre:
                fre_data += 1 #降水次数+1
                num_data += num #增加降水时数
                pre_array = np.append(pre_array,one_pre) #将单次降水量加入数组
            one_pre = 0
            numm = 0
            num = 0
    if fre_data == 0:
        mean_mun = 0
        precent = 0
        mean_data = 0
        max_pre = 0
        intensity = 0
    else:
        mean_mun = num_data/fre_data #平均每次的极端降水时长
        precent = dura_data/np.sum(array)#极端降水在总降水的贡献率
        max_pre = np.max(pre_array) #单次最大降水量
        mean_data = dura_data/fre_data #平均单次降水量
        intensity = dura_data/num_data #极端降水强度
    return dura_data,num_data,mean_num,precent,max_pre,mean_data,intensity #总极端降水量，时长...
                        
#根据极端降水事件定义，计算年降水事件特征---------------------------------------------------------------
def duration(array):
    row,col,band = array.shape
    dura_array=np.zeros((row,col))
    num_array=np.zeros((row,col))
    mean_num_array=np.zeros((row,col))
    precent_array=np.zeros((row,col))
    max_pre_array=np.zeros((row,col))
    mean_array=np.zeros((row,col))
    intensity_array=np.zeros((row,col))
    #创建空数组
    for i in range(row):
        for j in range(col):
            tempor=array[i,j,:].reshape(band)
            if tempor[0] == -999:
                mean_data = -999
            else:
                dura_data,num_data,mean_num,precent,max_pre,mean_data,intensity = DefineEx(tempor,band)
            dura_array[i,j]=dura_data
            num_array[i,j]=num_data
            mean_num_array[i,j]=mean_num
            precent_array[i,j]=precent
            max_pre_array[i,j]=max_pre
            mean_array=[i,j]=mean_data
            intensity_array[i,j]=intensity
    return dura_array,num_array,mean_num_array,precent_array,max_pre_array,mean_array,intensity_array

yearly_data=[]
#设置循环--------------------------------------------------------------------------------------------
for i in range(2001,2022):
    print(i)
    for n in range(1,13):
        #读取文件,将一年的数据存为一个数组
        newpath = os.path.expanduser('~/gee_data')
        newpath = os.path.join(newpath,str(i)+"/Hour_data"+str(i)+"_"+str(n)+".npy")
        b= np.load(file = newpath)
        #获取小时数据
        if n == 1:
            yearly_data = b.copy()
        else:
            yearly_data = np.concatenate((yearly_data,b),axis=2)
    #调用函数，计算特征
    if i == 2001:
        dura_array,num_array,mean_num_array,precent_array,max_pre_array,mean_array,intensity_array = duration(yearly_data)
        row,col=dura_array.shape
        dura_array=dura_array.reshape(row,col,1)
        num_array=num_array.reshape(row,col,1)
        mean_num_array=mean_num_array,.reshape(row,col,1)
        precent_array=precent_array.reshape(row,col,1)
        max_pre_array=max_pre_array.reshape(row,col,1)
        mean_array=mean_array.reshape(row,col,1)
        intensity_array=intensity_array.reshape(row,col,1)
        
    else:
        a,b,c,d,e,f,g = duration(yearly_data)
        dura_array = np.insert(dura_array,i-2001,a,axis=2)
        num_array = np.insert(num_array,i-2001,b,axis=2)
        mean_num_array = np.insert(mean_num_array,i-2001,c,axis=2)
        precent_array = np.insert(precent_array,i-2001,d,axis=2)
        max_pre_array = np.insert(max_pre_array,i-2001,e,axis=2)
        mean_array = np.insert(mean_array,i-2001,f,axis=2)
        intensity_array = np.insert(intensity_array,i-2001,g,axis=2)

#分别保存-----------------------------------------------------------------------------
newpath=os.path.expanduser('~/gee_data')
np.save(file=os.path.join(newpath, "dura_array.npy"), arr=dura_array)
np.save(file=os.path.join(newpath, "num_array.npy"), arr=num_array)
np.save(file=os.path.join(newpath, "mean_num_array.npy"), arr=mean_num_array)
np.save(file=os.path.join(newpath, "precent_array.npy"), arr=precent_array)
np.save(file=os.path.join(newpath, "max_pre_array.npy"), arr=max_pre_array)
np.save(file=os.path.join(newpath, "mean_array.npy"), arr=mean_array)
np.save(file=os.path.join(newpath, "intensity_array.npy"), arr=intensity_array)