GEE（python）雨天Gini指数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#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)