GEE（python）雨天Gini指数

一些碎碎念：最近一直忙着写论文，已经有大半年没有更新博文了。接下来几天会对过去半年写的代码做一些整理。希望文章能送审啊ballball了QAQ

雨天Gini指数是Gini指数的改进，用于分析降水在时间上的分布，也就是降水集中度。通过Gini指数(描述降水量在一年中的均匀分布情况)和年雨天数，使用Mann-Kendall检验和回归分析来评估两个指数的变化方向和速率。
来自于：Rajah, K., O’Leary, T., Turner, A., Petrakis, G., Leonard, M., & Westra, S. (2014). Changes to the temporal distribution of daily precipitation. GEOPHYSICAL RESEARCH LETTERS, 41(24), 8887-8894. http://doi.org/10.1002/2014GL062156

以下代码是利用GPM 2019-2021年(年份可更改)降水数据，将半小时时间尺度转为日尺度，计算雨天Gini指数，最终得每年的降水集中度数据，输出类型为numpy，可以通过gdal转为图像可视化。
（有关gdal的tif输入输出函数，可以参考我之前的博文python（GDAL）读取、输出tif数据 注：需提前获悉使用的降水数据的投影信息，可以通过在GEE里getInfo查询crs_transform；也可以先下载下来一张影像，用gdal读取投影信息。）

#Edited by Xinglu Cheng 2022.03.25
#读取半小时降水数据，转为日数据，并计算雨天GN指数
import ee
import os
import math
import threading
import numpy as np
import scipy.interpolate
import scipy.integrate
import matplotlib.pyplot as plt
from osgeo import gdal
from ipygee import Map
from IPython.display import Image
import sys

#加入研究区
HHH = ee.FeatureCollection("users/2210902126/GPM/HHH_Area")#该研究区为本人上传至GEE Assets的私人数据，需更改

#重采样-------------------------------------------------------------------------------------------
def repro(image):
    image=image.reproject(crs='EPSG:4326',scale=11132)
    return image

#获取影像的经纬度---------------------------------------------------------------------------------
def getImgData(image):
  image = image.addBands(ee.Image.pixelLonLat())
  data = image.reduceRegion(
      reducer = ee.Reducer.toList(),
      geometry = HHH,
      scale = 11132,
      maxPixels = 1e13
  )
  lat = np.array(ee.Array(data.get("latitude")).getInfo())
  lon = np.array(ee.Array(data.get("longitude")).getInfo())
  return lat, lon

#获取各种各样的信息-------------------------------------------------------------------------------
example=ee.ImageCollection("NASA/GPM_L3/IMERG_V06").filter(ee.Filter.date(ee.Date('2019-09-03').getRange('day'))).select('precipitationCal').sum().multiply(0.5).clip(HHH)#获取一个图像
example=repro(example)#重采样
print(example.getInfo())
row=geemap.ee_to_numpy(example,region=HHH,default_value=-999).squeeze().shape[0]
col=geemap.ee_to_numpy(example,region=HHH,default_value=-999).squeeze().shape[1]#获取行列号
lat,lon=getImgData(example)
lat_min=lat.max()
lon_min=lon.min()
print(lat_min,lon_min)


#将半小时数据变为日数据（用递归的方式）-------------------------------------------------------------
def toDaily(year, EndYear,DailyArray,bands):
    for i in range(1, 13):
        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):
            if i >= 10 and day >= 10: extend = ee.Date(str(year) + '-' + str(i) + '-' + str(d))
            if i >= 10 and day < 10: extend = ee.Date(str(year) + '-' + str(i) + '-0' + str(d))
            if i < 10 and day >= 10: extend = ee.Date(str(year) + '-0' + str(i) + '-' + str(d))
            if i < 10 and day < 10: extend = ee.Date(str(year) + '-0' + str(i) + '-0' + str(d))
            # 判断月份是否是两位数，两位数不需要加0
            image=ee.ImageCollection("NASA/GPM_L3/IMERG_V06").filter(ee.Filter.date(extend.getRange('day'))).select('precipitationCal').sum().multiply(0.5).clip(HHH).set({'system:time_start': extend})
            image=repro(image)#重采样
            
            if bands==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,bands,numArray,axis=2)#将数组存入
            bands+=1
            # 将半小时数据变为日数据,并加入日期属性,将数据加入数组
    year+=1
    if year > EndYear:
        return DailyArray,bands
    return toDaily(year, EndYear,DailyArray,bands)  # 递归

#定义GINI函数-------------------------------------------------------------------------------------
def gini(x, w=None):
    # The rest of the code requires numpy arrays.
    x = np.asarray(x)
    if w is not None:
        w = np.asarray(w)
        sorted_indices = np.argsort(x)
        sorted_x = x[sorted_indices]
        sorted_w = w[sorted_indices]
        # Force float dtype to avoid overflows
        cumw = np.cumsum(sorted_w, dtype=float)
        cumxw = np.cumsum(sorted_x * sorted_w, dtype=float)
        return (np.sum(cumxw[1:] * cumw[:-1] - cumxw[:-1] * cumw[1:]) / 
                (cumxw[-1] * cumw[-1]))
    else:
        sorted_x = np.sort(x)
        n = len(x)
        cumx = np.cumsum(sorted_x, dtype=float)
        # The above formula, with all weights equal to 1 simplifies to:
        return (n + 1 - 2 * np.sum(cumx) / cumx[-1]) / n

#去除无降水的数值---------------------------------------------------------------------------------
def pretreat(array):
    if (array < 0.1).any():
        b = np.where(array < 0.1)
        array = np.delete(array,b)
    return array

#应用函数-----------------------------------------------------------------------------------------

GData = np.array([])#创建一个空的数组存储结果
year = 0#计数
for i in range(2019,2022):
    bands=0#统计参与循环的图像数
    DailyArray=np.array([])#创建一个空的数组存储该年的日数据
    DailyArray,bands=toDaily(i,i,DailyArray,bands)#应用函数
    HArray=np.zeros((bands))#创建一个空的数组存储H(x)
    GArray=np.zeros((row,col))#用于存储该年的GINI值
    for n in range(row):
        for m in range(col):
            if DailyArray[n,m,1] == -999:
                Data = -999#剔除无效值
            else:
                tryArray = DailyArray[n,m,:].reshape(bands)
                tryArray = pretreat(tryArray)
                Data = gini(tryArray, w=None)
            GArray[n,m]=Data
    if year == 0:
        GData = GArray.copy().reshape(row,col,1)
    else:
        GData = np.insert(GData,year,GArray,axis=2)
    print(year+2001)
    year += 1  
print(GData.shape)