用python计算极端r95极端降水

徵言片羽

# 利用小时降水数据，统计西南五省市78个站点1981-2017年极端降水的频数和强度
# 西南五省市78个站点
import numpy as np
import math

#定义一个百分比函数，得到任意分位大于0小于2000的排位后的百分位数
#define a percentile function
def quantile_p(data:list, p):
    data.sort()
    ndata = []
    for i in data:
        if i>0 and i<2000:
            ndata.append(i)
    q = np.percentile(ndata,p)
    #print(ndata)
    return q

st1 = [56144,56146,56172,56178,56182,56196,56247,56251,56287,56374,\
     56385,56444,56459,56462,56475,56479,56485,56492,56533,56548,\
     56571,56586,56651,56664,56671,56684,56691,56739,56748,56751,\
     56763,56768,56778,56786,56793,56856,56872,56875,56886,57206,\
     57237,57306,57313,57328,57348,57405,57426,57432,57516,57606,\
     57614,57633,57647,57707,57718,57722,57729,57731,57741,57803,\
     57806,57816,57825,57832,57839,57902,57906,57910,57916,57922,\
     57932,57947,59007,59021,59023,59027,59046]
st = [56385]
dir  = 'D:/data/hourly_precip/1981-2017/PRE_hourly_'
endf = '.txt'
#prc = [] #设置空列表，存放降水值
fout = open('./PRE_hourly_frequency_itensity.txt','w')
for i in range(0,len(st)):
    filename = dir+str(st)+endf  #拼接构造文件绝对路径
    #print(filename)
    '''下面两行代码将多维数组转换为一维列表
     mulArrays = [[1,2,3],[4,5,6],[7,8,9]]
     list(np.array(mulArrays).flatten())
    '''
    #1.读取站点和经纬度信息
    fp = open(filename,'r')     #打开小时降水数据
    lines = fp.readlines()
    l1 = lines[0]
    stid = int(l1.split( )[0])      #站点信息
    lat  = int(l1.split( )[1])     #纬度，单位：度，分
    latitude = math.modf(lat/100)[1] + math.modf(lat/100)[0]*100/60  #纬度，单位：度

    #print('stid:',stid,'latitude:',latitude)
    lon  = int(l1.split( )[2])      #经度，单位：度，分
    longitude= math.modf(lon/100)[1] + math.modf(lon/100)[0]*100/60  #经度，单位：度


    #2.计算r95,r99等
    percent = 95  #极端降水指数，1-100任意值，但
    a = np.loadtxt(filename)    #将文本数据读入矩阵
    b = a[:,7:31]               #小时降水数据，7-31列，m*24维
    prc = list(np.array(b).flatten())  #矩阵转换为一维列表
    prc.sort()                  #列表排序
    rp = quantile_p(prc,percent)    #引用百分比函数计算当前站点降水r95阈值
    rp_prc = rp/10  #r95的真实阈值，单位mm

    #3.统计r95降水的频率frequency，用一维矩阵prc和r95的大小来判断
    n = 0
    for p in prc:
        if p>=rp:
            n+=1
    print('stid:',stid,'latitude:%.2f'%latitude,'longitude:%.2f'%longitude,'r95p:%smm'%rp,'frequency:',n)

徵言片羽

最近在做极端降水，需要计算r95，无奈没有现成现成脚本，就自己写了一个，供大家参考。重点就在于使用math自带的percentile函数，自己构造一个函数，函数对列表进行排序，选择其中大于0小于2000的数据，输出为新列表，再用math.percentile函数计算任意分位数，可以计算r95，r99等。不足和错误之处请斧正。

feiji158

厉害了，学习学习

又又又加零

顶顶 最近用NCL也做了极端指数 学习下

Polariscoco

学习啦

南瓜是我

码住学习

气象小白学编程

你好，我想问下为什么要选择其中大于0小于2000的数据呢。