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 = 31if i == 4 or i == 6 or i == 9 or i == 11: day = 30if i == 2 and year % 4 != 0 or i == 2 and year % 4 == 0 and year % 400 != 0: day = 28if 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))# 判断月份是否是两位数，两位数不需要加0image=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+=1if year > EndYear:return DailyArray,bandsreturn 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 overflowscumw = 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]=Dataif 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)

GEE（python）雨天Gini指数相关推荐

python 身体BMI指数判断
3.30 python身体BMI指数判断代码: height,weight=eval(input("请输入身高(米)体重(公斤):")) bim=weight/pow(heigh ...
Google Earth Engine（GEE）——关于NDVI_NDWI_NDBI_EVI指数波段运算公式函数的集成和优化（2）
上一篇文章写到了,函数合并处理的问题,链接如下: (17条消息) Google Earth Engine(GEE)--关于NDVI_NDWI_NDBI指数波段运算公式函数的集成和优化(1)_此星光明的 ...
【GEE python】基于geemap开发Google earth App
文章目录前言一.GEE python学习的"磨刀石" 二.geemap及其基本用法 1.geemap介绍 2. 比较常用的类 3.geemap的基本操作三.地图圈取,获取经纬 ...
Python float输出指数形式和小数形式切换（即科学计数法和完整数值切换）
Python float输出指数形式和小数形式切换(即科学计数法和完整数值切换) (本文数字全部瞎编,仅作举例之用,请勿相信. 另,本人Python初学者,写个文章作为笔记,也希望可以对他人有所帮助, ...
Python 二次指数平滑法预测
1.先上代码: import numpy as npdef secondaryExponentialSmoothingMethod(list,n_average,alpha,day): # 参数lis ...
Python 计算EMA(指数移动平均线)
总结使用递归和循环两种方法来完成 python环境下循环相比于递归更快,更适应极端样本情况递归 def _ema(arr,i=None):N = len(arr) α = 2/(N+1) #平滑指 ...
（87）--Python数据分析：指数密度函数与指数分布图
# 指数密度函数与指数分布图 lambd = 0.5 x = np.arange(0,15,0.1) y = lambd*np.exp(-lambd*x) plt.plot(x,y) plt.titl ...
python获得百度指数脚本[免费分享]
注意更新(2022-07-01日更新) 1. 估计是百度指数修改了爬虫策略,目前已更新为最新版本- 前言有时候大家需要知道一个关键词在互联网上的热度,想知道某个关键词的热度变化趋势.大家可能就是使 ...
python分析股票数据的项目_用Python分析股市指数
專欄 ❈本文作者:王勇,目前感兴趣项目商业分析.Python.机器学习.Kaggle.17年项目管理,通信业干了11年项目经理管合同交付,制造业干了6年项目管理:PMO,变革,生产转移,清算和资产处 ...
Google Earth Engine（GEE）——关于NDVI_NDWI_NDBI指数波段运算公式函数的集成和优化（1）
之前写过一篇文章是关于各种指数进行计算的文章,链接如下: (230条消息) 利用GEE(Google Earth Engine)在线处理NDVI.EVI.SAVI.NDMI等指数归一化教程!_此星光明 ...

GEE（python）雨天Gini指数

GEE（python）雨天Gini指数相关推荐

最新文章

热门文章