爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
12
返回列表 发新帖
楼主: cwllsss

[经验总结] 使用python进行一些气候相关统计量的计算

[复制链接]

新浪微博达人勋

 楼主| 发表于 2020-9-24 14:15:02 | 显示全部楼层
yya913 发表于 2020-9-24 08:30
用describe()就可以将mean,variance,百分位数等全部计算出来。但在处理气象数据的时候,比较麻烦的是缺测 ...

好呢谢谢你噢!会当心的!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-9-24 14:27:27 | 显示全部楼层
谢谢大家的很多关心、鼓励和帮助。还有yya913朋友的提醒!我知道python里面已经有很多数据处理的库以及函数可以非常方便地计算所需要的东西。但是因为我是在学习气候统计这个课程,所以自己就慢一点,想着尽量靠着自己的理解来把比较详细的过程都编写一遍,一是想靠着这个过程多熟悉熟悉python,二是为了加深对课程里面一些东西的理解。
真的很谢谢大家~~~本来以为就自己默默记录一下学习过程,没想到还有很多朋友会来关心~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-9-24 15:53:54 | 显示全部楼层
本帖最后由 cwllsss 于 2020-9-24 15:57 编辑

百分位数的计算
1、百分位数和百分位
百分位:将一组个数为n的数据值{x1,x2,x3,x4.......xn}从小到大排序。
                排序后为:{x(1),x(2),x(3),x(4).......x(n)}
则第p百分位,意为:在所有测量数据中,该测量值的累计频数达p%。
  • 累积频数(Cumulative Frequency),累积频数
    就是将各类别的频数逐级累加起来。通过累积频数,可以很容易
    看出某一类别(或数值)以下及某一类别(或数值) 以上的频数之和。

第p百分位所对应的测量值,称为百分位数。
  • 以身高为例,第5%百分位,表示有5%的人的身高小于此测量值,
    而有95%的人的身高大于此测量值。

2、几个常用的百分位数
  • 第25百分位数又称第一个四分位数(First Quartile),用Q1表示;
  • 第50百分位数又称第二个四分位数(Second Quartile),用Q2表示;
    第50百分位数就是中位数,与数据中心趋势有关。
  • 第75百分位数又称第三个四分位数(Third Quartile),用Q3表示。
  • 10%和90%分位数等与极值有关

3、百分位数的计算
    numpy.quantile
  • numpy.``quantile(a, q, axis=None,interpolation='linear' )
  • array_like
    Input array or object that can be converted to an array.
  • q
    array_like of float
    Quantile or sequence of quantiles to compute, which must be between 0 and 1 inclusive.
  • interpolation{‘linear’, ‘lower’, ‘higher’, ‘midpoint’, ‘nearest’}
    This optional parameter specifies the interpolation method to use when the desired quantile lies between two data points i < j:

    • linear: i + (j - i) * fraction, where fraction is the fractional part of the index surrounded by i and j.
    • lower: i.
    • higher: j.
    • nearest: i or j, whichever is nearest.
    • midpoint: (i + j) / 2.




方法一
第1步:以递增顺序排列原始数据(即从小到大排列)。
第2步:计算指数i=np%
第3步:
l)若 i 不是整数,将 i 向上取整。大于i的毗邻整数即为第p百分位数的位置。
2) 若i是整数,则第p百分位数是第i项与第(i+l)项数据的平均值。
import numpy as np
import pandas as pd
import math
'''
#整数、小数分离
math.modf(num)
#向上取整
math.ceil(num)
'''
def quantile(arr,p):
    p_quantile = 0
    #从小到大排序
    arr_sort = np.sort(arr)
    n = np.size(arr_sort)
    #计算第p百分位的项数
    i = n*(p/100.0)
    print(i)
    if(math.modf(i)[1] == 0):
        p_quantile = (arr_sort[i-1]+arr_sort)/2.0
    else:
        p_quantile = arr_sort[math.ceil(i)]
    return p_quantile
               
#生成随机数
arr1 = np.random.randint(1,100,size=111)
print(arr1)
print(np.sort(arr1))
#第50百分位数
p_quantile1 = quantile(arr1,50)
#中位数
medianarr1 = np.median(arr1)
print(p_quantile1)
print(medianarr1)
print(np.quantile(arr1,0.5))

方法二
第一步:将n个变量值从小到大排列,X(j)表示此数列中第j个数。
第二步:计算指数,设(n+1)P%=j+g,j为整数部分,g为小数部分。
第三步:1)当g=0时:P百分位数=X(j);
2)当g≠0时:P百分位数=gX(j+1)+(1-g)X(j)=X(j)+g*[X(j+1)-X(j)]。

import numpy as np
import pandas as pd
import math
'''
#整数、小数分离
math.modf(num)
#向上取整
math.ceil(num)
'''
def quantile(arr,p):
    p_quantile = 0
    #从小到大排序
    arr_sort = np.sort(arr)
    n = np.size(arr_sort)
    #计算指数,设(n+1)P%=i_int+i_float,i_int为整数部分,i_float为小数部分。
    i = (n+1)*(p/100.0)
    #分离出i的小数、整数部分
    i_float = math.modf(i)[0]
    i_int = int(math.modf(i)[1])
    print(i,i_float,i_int)
   
    #当i_float=0时
    if(i_float == 0):
        p_quantile = arr_sort[i_int]
        
    #当i_float≠0时
    else:
        p_quantile = arr_sort[i_int-1]+i_float*(arr_sort[i_int] - arr_sort[i_int-1])
    return p_quantile
        
#生成随机数
arr1 = np.random.randint(1,100,size=78)
print(arr1)
print(np.sort(arr1))
#第50百分位数
p_quantile1 = quantile(arr1,50)
#中位数
medianarr1 = np.median(arr1)
print(p_quantile1)
print(medianarr1)
print(np.quantile(arr1,0.5))




密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-9-26 13:07:39 | 显示全部楼层
Keep going, good job
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-9-27 15:46:13 | 显示全部楼层
加油,我是刚刚开始学习Python语言的气象小学生,向您学习,请不吝赐教!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-9-27 20:34:30 | 显示全部楼层
VV09 发表于 2020-9-27 15:46
加油,我是刚刚开始学习Python语言的气象小学生,向您学习,请不吝赐教!

一起学习!!一起加油!!!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2020-10-2 17:04:57 | 显示全部楼层
本帖最后由 cwllsss 于 2020-10-3 11:23 编辑

距平-Anomaly

定义式

x = xi- x &#773;


    import xarray as xr
    import numpy as np
    import pandas as pd
    import os
    os.environ['PROJ_LIB'] = r'D:\anaconda3\pkgs\share\proj'
    from sklearn.feature_selection import f_regression
    import cartopy.feature as cfeat
    import matplotlib.pyplot as plt
    #import matplotlib;matplotlib.use('TkAgg')
    from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
    from mpl_toolkits.basemap import Basemap
   
    from IPython.core.pylabtools import figsize # import figsize
    #path = os.getcwd()#获取当前路径
    #print(path)
    #os.chdir('/disk1/cwl/data')
   
    #读取数据
    foo = xr.open_dataset('D:\postgraduate study\data\sst.mon.mean.2.5.nc')
    #print(foo)
   
    #取出全时段sst数据
    sst_all = foo['sst']
   
    #取出1979-2018时段的数据
    sst = sst_all.loc[:'2018-12-01',:,:]
   
    #print(sst)
   
    #求1979-2018四十年平均
    sst_40ymean = np.mean(sst,0)
   
    #求2009-2018十年平均
    sst_10y = sst.loc['2009-01-01':,:,:]
    sst_10ymean = np.mean(sst_10y,0)
   
    #2009-2018相对于1979-2018四十年平均
   
    anomaly = sst_10ymean - sst_40ymean
   
    '''
    print(sst_40ymean)
    print(sst_10ymean)
    print(anomaly)
    #max and min
    print(np.max(anomaly))
    print(np.min(anomaly))
   
    sstamax = np.max(anomaly)
    sstamin = np.min(anomaly)
    '''
    lon_need = anomaly['longitude']
    lat_need = anomaly['latitude']
   
   
    #画图
    #区域设置( llcrnrlon 是起始经度, llcrnrlat 是起始纬度; urcrnrlon 是起始纬度; urcrnrlat 是终止纬度)
    map = Basemap(llcrnrlon=0, llcrnrlat=-70, urcrnrlon=180, urcrnrlat=90)
    #画海岸线和国界
    map.drawcoastlines()
    Lon,Lat= np.meshgrid(lon_need,lat_need)
    X,Y = map(Lon,Lat)
    lvls = np.arange(-1.6,1.7, 0.15)
    #画填色图
    C = map.contourf(Lon,Lat,anomaly,cmap = plt.cm.RdYlBu_r,levels = lvls,extend='both' )#                                            
    cbar = map.colorbar(C,'right',ticks = np.arange(-1.6,1.7, 0.15) )#extend='both' ,format = '%.1f'
   
    #设置经纬网格线
    parallels = np.linspace(-70,90,9)
    map.drawparallels(parallels,labels=[True,False,False,False])
    meridians = np.linspace(0,180,6)
    map.drawmeridians(meridians,labels=[False,False,False,True])
    #设置标题
    plt.title(r'SST Anomaly(K)',fontsize=12)
    #设置图片大小和分辨率
    plt.rcParams['savefig.dpi'] = 300 #图片像素
    plt.rcParams['figure.dpi'] = 300 #分辨率
    #保存图片
    plt.savefig('D:\postgraduate study\data\ssta1.png')
———————————————————————————————————————————————————————
1.png
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-1-7 09:54:27 | 显示全部楼层
{:eb502:}
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-1-18 17:20:12 | 显示全部楼层
{:eb502:}{:eb502:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表