1.引用库2.定义函数4.趋势分析及MK检验5.结果图像保存至本地
这段时间有关GPM降水分析的一些代码总结：
1.引用库#Edited by Xinglu Cheng 2022.01.06
#读取半小时降水数据变为日数据,获取超过90%阈值的总降水量
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")
2.定义函数
#归一化处理---------------------------------------------------------------------------------------------------------------------------
def UnitToOne(image):
minMax = image.reduceRegion(
reducer = ee.Reducer.minMax(),
geometry = HHH,
scale = 11132,
maxPixels = 10e9
)
name = ee.String(image.bandNames().get(0))
return image.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))))
#重采样-------------------------------------------------------------------------------------------
def repro(image):
image=image.reproject(crs='EPSG:4326',scale=11132)
return image
#半小时数据转为日数据函数--------------------------------------------------------------------------
DailyPre=ee.List([])
def toDaily(year, EndYear, DailyPre):
if year > EndYear:
return DailyPre
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
DailyPre= DailyPre.add(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}))
# 将半小时数据变为日数据,并加入日期属性
year+=1
return toDaily(year, EndYear, DailyPre) # 递归
#获取影像的经纬度---------------------------------------------------------------------------------
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)#输出左上角经纬度
3.调用函数以及reducer进行分析
方法1的结果可以直接用geemap进行显示；
方法2由于数据量太大，转为numpy处理，用gdal转为图像存至本地。
#方法1：对每个像元进行筛选，统计超过百分比阈值的日数------------------------------------------------
# year_begin = 2001
# year_end = 2002
# DailyPre=toDaily(year_begin,year_end,DailyPre)#应用函数
# DailyCollection=ee.ImageCollection.fromImages(DailyPre)#将list变为图像集
# def subtract(image):
# image=image.subtract(Percentile)#减去阈值
# image=image.mask(image.gt(0))#把小于0的部分变为无效值
# return image
# Percentile = DailyCollection.reduce(ee.Reducer.percentile([90]))#计算区间百分比内的平均值
# GreaterSubtract = DailyCollection.map(subtract)
# CountDaily= GreaterSubtract.reduce(ee.Reducer.count())#计算每个像元非零值的个数，即大于阈值的个数
#方法2：计算20年的极端降水的总降水量变化趋势--------------------------------------------------------
#定义函数：将阈值以下的部分掩膜
def calculate(year,DailyPre):
DailyPre=toDaily(year,year,DailyPre)#应用函数
DailyCollection=ee.ImageCollection.fromImages(DailyPre)
Percentile = DailyCollection.reduce(ee.Reducer.percentile([90]))#计算区间百分比值
def subtract(image):
imageb=image.subtract(Percentile)#减去阈值
image=image.mask(imageb.gt(0))#把小于0的部分变为无效值
return image
DailySubtract = DailyCollection.map(subtract).sum().set("date",ee.Number(i))#提取极端降水事件，计算极端降水的总降水量,并添加年份属性
return DailySubtract
year_begin = 2001
year_end = 2022
bands = 0
for i in range(year_begin,year_end):
DailySubtract = calculate(i,DailyPre)
DailySubtract = repro(DailySubtract)
if bands==0:
numArray=geemap.ee_to_numpy(DailySubtract,region=HHH,default_value=-999)
DailyArray=numArray.copy()
else:
numArray=geemap.ee_to_numpy(DailySubtract,region=HHH,default_value=-999).squeeze()#图像转数组
DailyArray=np.insert(DailyArray,bands,numArray,axis=2)#将数组存入
bands += 1
print(i)
4.趋势分析及MK检验
使用方法2的结果数组（每年超过90%百分值的总降水量）进行分析
#线性倾向-----------------------------------------------------------------------------------------
# from sklearn.linear_model import LinearRegression
# #使用线性回归模型进行建模
# y = np.arange(2001,2022,1).reshape(-1, 1)#因变量
# LinearData = np.array([])#创建一个空的数组存储数据
# for i in range(row):
# for j in range(col):
# x = DailyArray[i,j,:].reshape(bands).reshape(-1, 1)
# if x[0] == -999:
# Data = -999
# else:
# LinearModel = LinearRegression()
# LinearModel.fit(x,y)#使用自变量x和因变量y进行训练模型
# Data = float(LinearModel.coef_[0])
# LinearData = np.append(LinearData,Data)
# LinearData=LinearData.reshape(row,col)#变为二维数组
# print(LinearData)
#标准差，中位数 验证波动性-------------------------------------------------------------------------
# LinearData = np.array([])#创建一个空的数组存储数据
# for i in range(row):
# for j in range(col):
# x = DailyArray[i,j,:].reshape(bands)
# if x[0] == -999:
# Data = -999
# else:
# Data = np.std(x)#求标准差
# # Data = np.median(x)#求中位数
# LinearData = np.append(LinearData,Data)
# LinearData=LinearData.reshape(row,col)#变为二维数组
# print(LinearData)
#sen_MK趋势分析-----------------------------------------------------------------------------------
import pymannkendall as mk
LinearData = np.array([])#创建一个空的数组存储数据
for i in range(row):
for j in range(col):
x = DailyArray[i,j,:].reshape(bands)
if x[0] == -999:
Data = -999
else:
trend, h, p, z, Tau, s, var_s, slope, intercept = mk.original_test(x)
#trend:增加减少或没有趋势（1.，-1，0）
#h:趋势是否存在，True为存在，False为不存在
#p:显著性检验的p值
#z:normalized test statistics 正态分布检验统计
#Tau:Kendall Tau 相关系数
#s:Mann-Kendal's score
#var_s:Variance S
#slope:senslope
if trend=="decreasing":
trend_value=-1
elif trend=="increasing":
trend_value=1
else:
trend_value=0
Data=slope#设置输出的要素
LinearData = np.append(LinearData,Data)
LinearData=LinearData.reshape(row,col)#变为二维数组
print(LinearData)
#MK突变分析---------------------------------------------------------------------------------------
# from scipy import stats
# def sk(data):
# n=len(data)
# Sk = [0]
# UFk = [0]
# s = 0
# E = [0]
# Var = [0]
# for i in range(1,n):
# for j in range(i):
# if data[i] > data[j]:
# s = s+1
# else:
# s = s+0
# Sk.append(s)
# E.append((i+1)*(i+2)/4 ) # Sk[i]的均值
# Var.append((i+1)*i*(2*(i+1)+5)/72 ) # Sk[i]的方差
# UFk.append((Sk[i]-E[i])/np.sqrt(Var[i]))
# UFk=np.array(UFk)
# return UFk
# #a为置信度
# def MK(data,a):
# ufk=sk(data) #顺序列
# ubk1=sk(data[::-1]) #逆序列
# ubk=-ubk1[::-1] #逆转逆序列
# #输出突变点的位置
# p=[]
# u=ufk-ubk
# for i in range(1,len(ufk)):
# if u[i-1]*u[i]<0:
# p.append(i)
# if p:
# return p
# else:
# p = 0
# return p
# #画图
# # conf_intveral = stats.norm.interval(a, loc=0, scale=1) #获取置信区间
# # plt.figure(figsize=(10,5))
# # plt.plot(range(len(data)),ufk,label = 'UFk',color = 'r')
# # plt.plot(range(len(data)),ubk,label = 'UBk',color = 'b')
# # plt.ylabel('UFk-UBk',fontsize=25)
# # x_lim = plt.xlim()
# # plt.ylim([-6,7])
# # plt.plot(x_lim,[conf_intveral[0],conf_intveral[0]],'m--',color='r',label='95%显著区间')
# # plt.plot(x_lim,[conf_intveral[1],conf_intveral[1]],'m--',color='r')
# # plt.axhline(0,ls="--",c="k")
# # plt.legend(loc='upper center',frameon=False,ncol=3,fontsize=20) # 图例
# # plt.xticks(fontsize=25)
# # plt.yticks(fontsize=25)
# # plt.show()
# #输入数据和置信度即可
# LinearData = np.array([])#创建一个空的数组存储数据
# for i in range(row):
# for j in range(col):
# x = DailyArray[i,j,:].reshape(bands)
# if x[0] == -999:
# Data = -999
# else:
# Data=MK(x,0.95)#设置输出的要素
# if Data!=0:
# Data=len(Data)
# LinearData = np.append(LinearData,Data)
# LinearData=LinearData.reshape(row,col)#变为二维数组
# print(LinearData)
5.结果图像保存至本地
#输出图像
def write_tif(newpath,im_data,im_Geotrans,im_proj,datatype):
if len(im_data.shape)==3:
im_width, im_height, im_bands = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
diver = gdal.GetDriverByName('GTiff')
new_dataset = diver.Create(newpath, im_width, im_height, im_bands, datatype)
new_dataset.SetGeoTransform(im_Geotrans)
new_dataset.SetProjection(im_proj)
if im_bands == 1:
new_dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
new_dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del new_dataset
#设置投影和仿射变换
im_Geotrans=(lon_min, 0.10000045742818513, 0.0, lat_min, 0.0, -0.10000045742818513)#投影参数随图像修改
im_proj='GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
#调用输出函数
newpath=os.path.expanduser('~')
newpath = os.path.join(newpath, '_name_.tif')#输入名称
write_tif(newpath,LinearData,im_Geotrans,im_proj,gdal.GDT_Float32)
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）