NCL统计每个站气温年平均大于35°日数

用NCL统计了每个站气温大于35°的日数再处以年数，就可以得到每个站年平均发生EHE的次数了。代码如下：

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
    
    FillValue=-9999.0
    fin=addfile("Tmax.nc","r")
    T0=fin->Tmax                ;把Tmax写进T0里去，但是T0还需要除10才是真实的最高气温
    Tmax=T0*0.1                 ;至此Tmax是真正的每站每日最高气温了
    size=dimsizes(Tmax)
    Tmax=where(ismissing(Tmax),0,Tmax);我把缺测值都换成0了，这样在统计高温的时候就不会报错了
    ntime=size(0)     ;时间序列数
    nst=size(1)       ;站点数    
    delete(size)
    ; 开始计算每个站点最高温度超过35的日子，count用来计数
    count=new(nst,"float",0)
    do i=0,nst-1
        count(i)=num(Tmax(:,i).gt.35.0);数出所有年份发生极端高温的日数
    end do
    count=count/56    ;至此年平均EHE日数已经统计完毕

    file_path="/mnt/d/hong/documents/project/temp/SURF_CLI_CHN_MUL_DAY-TEM-12001-201612.txt"
    nrow=numAsciiRow(file_path)
    ncol=numAsciiCol(file_path)
    data=asciiread(file_path, (/nrow,ncol/), "float")

    lat_a=round(data(:,1)*0.01,0)+(data(:,1)*0.01-round(data(:,1)*0.01,3))/60*100
    lon_a=round(data(:,2)*0.01,0)+(data(:,2)*0.01-round(data(:,2)*0.01,3))/60*100
    lat=new(nst,"float")
    lon=new(nst,"float")
    do i=0,nst-1
        lon(i)=lon_a(31*i)  ;读入经度
        lat(i)=lat_a(31*i)  ;读入纬度
    end do

    R=count                ;要素水平就是年平均EHE日数
    arr=(/10,20,30,40,50/)
    colors = (/10,38,56,66,75,94/);5个水平需要5个颜色来代表
    num_markers=dimsizes(arr)+1;
    lat_new = new((/num_markers,dimsizes(R)/),float,-999);
    lon_new = new((/num_markers,dimsizes(R)/),float,-999);
    labels = new(dimsizes(arr)+1,string);最后出现在图下方的标签
    do i =0,num_markers-1
        if(i.eq.0);第一个水平即低于0的水平
            indexes=ind(R.lt.arr(0))
            labels(i)="R<"+arr(0)
        end if
        if(i.eq.(num_markers-1))then;最大的一个水平即为大于26的
            indexes=ind(R.ge.max(arr))
            labels(i)="R>="+max(arr)
        end if
        if(i.gt.0.and.i.lt.(num_markers-1))then;中间的水平
            indexes=ind(R.ge.arr(i-1).and.R.lt.arr(i))
            labels(i)=arr(i-1)+"<=R<"+arr(i)
        end if
        if(.not.any(ismissing(indexes)))then;如果index里有数，而不全是-999，那么把lat、lon_new的前N个（在这一水平里有N个点）换成这N个点的经纬度。
            npts_range=dimsizes(indexes)
            lat_new(i,0:(npts_range-1))=lat(indexes)
            lon_new(i,0:(npts_range-1))=lon(indexes)
        end if
        delete(indexes)
    end do
    wks=gsn_open_wks("png", "Number of EHE");文件名及格式
    gsn_define_colormap(wks, "WhViBlGrYeOrRe")
    mpres=True
    mpres@gsnFrame=False
    mpres@pmTickMarkDisplayMode="Always"
    mpres@mpMinLatF=15.0
    mpres@mpMaxLatF=55.0
    mpres@mpMinLonF=72
    mpres@mpMaxLonF=135.0
    mpres@tiMainString="Number of EHE";图形标题
    mpres@mpDataBaseVersion = "Ncarg4_1";这一步和下一步必须要，否则加载中国地图的时候会出错（找不到地图库）
    mpres@mpDataSetName="Earth..4";这步加上步再加下面那个China:state和Taiwan就可以画出中国轮廓边界了
    mpres@mpOutlineOn=True
    mpres@mpOutlineSpecifiers=(/"China:states","Taiwan"/);在这个地图库里我们绘制中国和台湾的区域边界，China:state里居然没有台湾！不能忍，要包括进来！
    mpres@mpGeophysicalLineThicknessF=2.0 ;这两行是为了加粗边界和国界线
    mpres@mpNationalLineThicknessF=2.0
    mpres@mpDataBaseVersion = "Ncarg4_1"
    mpres@mpAreaMaskingOn = True   ;使能填充覆盖
    mpres@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/)
    map=gsn_csm_map(wks,mpres)
    gsres=True                  ;添加标识别
    gsres@gsMarkerIndex=16      ;用实心圆形式
    txres=True                  ;gsn_text_ndc中最后参数
    txres@txFontHeightF=0.015   ;设置字体大小

    ;图例和文本的位置（单位坐标系内的坐标位置）
    xleg=(/0.07,0.07,0.32,0.32,0.57,0.57,0.82,0.82/)
    yleg=(/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/)
    xtxt=(/0.16,0.16,0.41,0.41,0.66,0.66,0.91,0.91/)
    ytxt=(/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/)

    ;画几条线分割北方地区
    lnres=True
    lnres@gsLineColor="black"
    lnres@gsLineThicknessF="5.0"
    lnres@gsLineDashPattern=0;实线
    str=unique_string("string")
    llon=(/80.3,119.0/)
    llat=(/35.0,35.0/)
    map@$str$=gsn_add_polyline(wks,map,llon,llat,lnres)
    draw(map);先画map和map上的一条线（我设置的35°N）后面才是画图例和标识点

    ;依次绘制各个数值区间范围的标识点及其图例
    do i=0,num_markers-1
        if (.not.ismissing(lat_new(i,0))) then
            gsres@gsMarkerColor=colors(i)
            ;gsres@gsMarkerThicknessF=0.5*i+0.5
            gsres@gsMarkerSizeF=5.0+1.0*i
            gsn_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres)
        
        ;绘制图形下方的标识和文本以作为图例
        gsn_polymarker_ndc(wks, xleg(i), yleg(i), gsres)
        gsn_text_ndc(wks,labels(i), xtxt(i), ytxt(i), txres)
        end if
    end do
    frame(wks);翻页

end