请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 10417|回复: 4

[经验总结] ncl空间分布图代码

[复制链接]

新浪微博达人勋

发表于 2019-6-5 16:16:12 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 veroniaca 于 2019-6-5 16:22 编辑

如果想要画出类似下面的图片,您可以好好看看如下代码:(黄色背景都是数据的处理,可以直接不看。绿色字是无关紧要的口令,可直接省略。其他均是绘图的设置)这个代码出来的图片,就是我下面贴出来的这张。
T2.PNG ·
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

;wrf模式资料的数据
          filelist=systemfunc("ls /mnt/d/wrf*")     ;读取多个文件名用filelist
          f5=addfiles(filelist,"r")                          ;读取数据,r表示只能读取
          temdata=wrf_user_getvar(f5,"T2",-1)    ;读取数据
          wrfLon=wrf_user_getvar(f5,"XLONG",1) ;读取数据
          wrfLat=wrf_user_getvar(f5,"XLAT",1)     ;读取数据
          lon5=fspan(99.75,125.25,35)                ;创建格点,从99.75到125.25之间,生成35个格点
          lat5=fspan(0,25.5,35)                           ;创建格点,从0到25之间,生成35个格点
          result=rcm2rgrid_Wrap(wrfLat,wrfLon,temdata,lat5,lon5,1) ;这里进行了差值运算,rcm2rgrid_Wrap是反距离权重函数,将原本的格点数据插值到创建的格点上
          ave5=dim_avg_n(result,0)                     ;对数据进行平均
          aave5=ave5-273.15
        ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ;设置经度属性(这里要设置属性,是因为插值后的lon5,lat5只是一个二维格点,没有属性)
        lon5!0                        ="lon"                            ;给lon5设置的属性是经度属性
        lon5@long_name           ="lon"
        lon5@units                ="degrees_east"             ;将北纬设置为正值,因此南纬为负值,比如50°E,表示-50
        lon5&lon                        =lon5
        ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ;设置纬度属性
        lat5!0                        ="lat"
        lat5@long_name            ="lat"
        lat5@units                        ="degrees_north"
        lat5&lat                        =lat5
        ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ;一维序列插值到二维网格,然后变量第一维设为纬度,第二维设为经度
        ;这是ncl画图默认一维纬度,二维经度
        R5                                =aave5                             ;这里是给要画图的变量设置属性
        R5!0                                ="lat"                               ;在ncl里面,从0开始数,因此这表示第一维
        R5!1                                ="lon"
        R5&lat                        =lat5
        R5&lon                        =lon5
;NCEP资料的数据
        filename1="/mnt/f/dasixia/data/ncep/air/air.2m.gauss.2015.nc"  ;读取文件名
        f1=addfile(filename1,"r")                ;打开文件,r表示只读
        temdata1=f1->air                          ;读取数据,将f1文件里的air变量,命名为temdata1
        lon1=f1->lon                                ;读取数据,将f1文件里的lon变量,命名为lon1
        lat1=f1->lat                                  ;读取数据,将f1文件里的lat变量,命名为lat1
        nlon1=dimsizes(lon1)                     ;dimsizes指的是数组的个数,这里用nlon1表示文件里lon1的个数
        nlat1=dimsizes(lat1)
        dim1=dimsizes(temdata1)                          ;读取要素的个数
        ave1=new((/dim1(1),dim1(2)/),"float")        ;new表示创建数组,其中确定的是数组的维数和种类(float指的是浮点数)。dim1(1)表示dim1的第二维,即经度或者纬度。这里表示生成的是二维数组
        ave1=0                                                     ;因为ave此时是垃圾数,将其所有值设为0
        copy_VarCoords(temdata1(1,:,:),ave1)         ;将temdata这个数组复制给ave这个数组,并且此时ave所具有的维数特征和temadata一样
        do i=273,300                                             ;这是时间。从nc数据可知,有一年的数据。1对应的是xx年1月1日。这里273对应的是2015年10月1日,这里只取了28天的数据
;             if(i.eq.197.or.i.eq.200.or.i.eq.210)then   ;选取缺测值,若是这几个天数,则忽略不计
;                     continue
;             end if

             ave1=ave1+temdata1(i,:,:)                    ;赋值,这个循环是将28天的数据相加
        end do
        ave1=ave1/28-273.15                                 ;除于28是因为前面数据是28天相加
;EC资料的数据
        filename2="/mnt/f/dasixia/data/ec/grib2.nc"     ;读取文件名,这里下载ERA-Interim再分析资料的2015年10月的每6小时数据
        f2=addfile(filename2,"r")                                 ;打开文件,r表示只读
        temdata2=f2->t2m                                         ;读取数据,将f2文件里的t2m变量,命名为temdata2
        add=temdata2@add_offset                              ;这个语句是ERA-interim数据独有的;add_offset是数据里的类似起到调整作用的数值
        scale=temdata2@scale_factor                          ;这个语句是ERA-interim数据独有的;scale_factor是数据里的类似起到调整作用的数值
        lon2=f2->longitude                                         ;读取数据,将f2文件里的longitude变量,命名为lon2
        lat2=f2->latitude                                            ;读取数据,将f1文件里的latitude变量,命名为lat2
        nlon2=dimsizes(lon2)                                      ;dimsizes求的是数组的个数,这里用nlon2表示文件里lon2的个数
        nlat2=dimsizes(lat2)
        dim2=dimsizes(temdata2)
        dat2=new((/dim2(0),dim2(1),dim2(2)/),"double")  ;new表示创建数组,其中确定的是数组的维数和种类(ERA-Interim数据的变量的记录是以double的形式记录)。这里表示生成的是三维数组
        dat2=temdata2*scale+add                                  ;赋值,这里是将所求的温度进行了处理,这是ERA-interim数据独有的且是必须的处理方式。主要是因为其有其他数据都没有的调整值
        ave2=new((/dim2(1),dim2(2)/),"double")              ;new表示创建数组,其中确定的是数组的维数和种类(ERA-Interim数据的变量的记录是以double的形式记录)。这里表示生成的是二维数组
        ave2=0                                                              ;因为ave2此时是垃圾数,将其所有值设为0
        aave2=new((/dim2(1),dim2(2)/),"double")
         aave2=0
        j=0
         copy_VarCoords(temdata2(1,:,:),aave2)                ;将temdata2这个数组复制给aave2这个数组,并且此时aave2所具有的维数特征和temadata2一样
        do i=0,111                                                          ;112/4=28天,因为数据是从2015年10月开始,一天有四个时次的资料,这里选用28天,因此是112个数据,但是由于ncl是从0开始数,因为是0-111
                ave2=ave2+dat2(i,:,:)                                 ;赋值,将dat2这个数组(即每个时间段的温度)都赋值到ave2
                j=j+1                                   
                if(j.eq.4)then                                               ;若j=4
                        aave2=ave2/4+aave2                           ;则将ave2除于4,这是为了求出一天的平均,再加上aave2,这是为了将28天的数据相加
                        ave2=0                          
                        j=0
                end if
        end do
         aave2=aave2/28-273.15                                      ;前面aave2是28天相加,因此需要除以28
;JRA-55资料的数据
         diri ="/mnt/f/dasixia/data/jra/anl_surf125.201510"
         files3 = systemfunc("ls "+ diri + "*.grb")
         grb =addfiles(files3,"r")
         temp3= new((/112,145,288/),float)
         temp3= 0
         do  i = 0,111                  ;112/4=28天的数据录入
         temp31=grb->TMP_GDS0_HTGL
         temp3(i,:,:)=temp31
         end do
         lon3 =grb[1]->g0_lon_1
         lat3 =grb[1]->g0_lat_0
         nlat3=dimsizes(lat3)        
         nlon3=dimsizes(lon3)                        
         dim3=dimsizes(temp3)        
         ave3=new((/dim3(1),dim3(2)/),"float")               
         ave3=0                                                                        
         copy_VarCoords(temp3(1,:,:),ave3)               
         do j=0,111        
                 ave3=ave3+temp3(j,:,:)
         end do
         ave3=ave3/112-273.15
; ;NOAA资料的数据
        file_path                        ="/mnt/f/dasixia/data/noaa/air/finalRes.txt"
        data                                =readAsciiTable(file_path,33,"float",1)
        lon                                 =data(:,0)
        lat                                =data(:,1)
        temData                                 =data(:,2:)
        temData@_FillValue                =32766
        temAve                                =new(1118,"float")
        temData=(temData-32)/1.8
        do i = 0,1117
                temAve(i)=avg(temData(i,:))
        end do

        lonNew = new(35,"float")
        latNew = new(35,"float")
        do j = 0,34,1
                lonNew(j)=j*0.75+99.75
        end do
        do k = 0,34,1
                latNew(k)= k*0.75
        end do
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ;设置经度属性
        lonNew!0                        ="lon"
        lonNew@long_name                ="lon"
        lonNew@units                        ="degrees_east"
        lonNew&lon                        =lonNew
        ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ;设置纬度属性
        latNew!0                        ="lat"
        latNew@long_name                ="lat"
        latNew@units                        ="degrees_north"
        latNew&lat                        =latNew
        ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        ;一维序列插值到二维网格,然后变量第一维设为纬度,第二维设为经度
        ;这是ncl画图默认一维纬度,二维经度
        R                                =cssgrid(lat,lon,temAve,latNew,lonNew)
        R!0                                ="lat"
        R!1                                ="lon"
        R&lat                                =latNew
        R&lon                                =lonNew


;-- define the workstation
        cmap                                        =read_colormap_file("NCV_blue_red")     ;引用的色板,详情可见http://www.ncl.ucar.edu/Document ... table_gallery.shtml
        wks_type                                   ="pdf"                  ;设置输出的图像文件格式为PDF,工作区设置就要按照pdf/ps的方式设置
;        wks_type@wkPaperWidth           =15                        ;设置画布的宽为15英寸
;        wks_type@wkPaperHeightF        =6                        ;设置画布的高为6英寸

        wks                                           =gsn_open_wks(wks_type,"/mnt/f/dasixia/data/T2")    ;图片输出的位置
        res                                            =True                   ;如果设置是正确
        res@lbLabelBarOn                      = False                 ;不显示色板,如果写True,则每一个图都会有对应的色板。
        res@gsnMaximize                       =True                   ;设置图像最大化,此时图像大小为画布大小减去边距
        res@gsnPaperOrientation             ="portrait"                ;设置图像的方向
        res@gsnPaperMargin                  =0.2                        ;设置图像到画布边缘的距离,即边距
        res@gsnAddCyclic                       =False                  ;如果设置是True,则循环点被加入数据;如果数据不是循环的,就设置为False就可以        
        ;res@gsnCenterString                 ="Temperature of the South China Sea in Oct. 2015"    ;设置图片的标题,标题会出现在每张图片的正中间
        res@mpDataSetName                 ="Earth..4"             ;中国地图包含在这个叫Earth..4的地图库里
        res@mpDataBaseVersion             ="MediumRes"       ;分辨率
        ;res@mpLimitMode                     ="LatLon"
        res@mpLimitMode                      ="Corners"           ;点对点的投影
        res@mpOutlineOn                      =True                  ;绘制地图区域轮廓
        res@mpLeftCornerLatF               =0                       ;设置图最南边的纬度
        res@mpRightCornerLatF             =25                     ;设置图最北边的经度
        res@mpLeftCornerLonF              =100                    ;设置图最西(左)边的经度
        res@mpRightCornerLonF            =125                    ;设置图最东(右)边的经度
        res@mpProjection                     ="Mercator"           ;投影方式,这里用的是墨卡托投影              
        res@mpDataBaseVersion            ="Ncarg4_1"
        ;res@cnRasterSmoothingOn       =True         
        res@cnFillOn                            =True                   ;是值能填充
        res@cnLinesOn                         =False                  ;是否画等值线
        res@cnLineLabelsOn                 =False                  ;不要等直线上的标签
        res@cnInfoLabelOn                   =False                  ;若设置为true,则表示局部最大值的区域放置标签。
        res@pmTickMarkDisplayMode    ="Conditional"          ;always效果一样
        res@tmXTOn                           =False                  
        res@tmYROn                           =False
        res@mpGridAndLimbOn            =False                  ;是否显示网格
        res@mpGridLatSpacingF           =2.5                    ;纬度网格的间距
        res@mpGridLonSpacingF          =2.5                    ;经度网格的间距
        res@mpGridLineThicknessF      =0.5                    ;线的粗细        
        res@cnFillPalette                    =cmap                       
        gsn_define_colormap(wks,"NCV_blue_red")           ;所引用的色板名称
        res@cnLevelSelectionMode="ExplicitLevels"           ;是否自己设置等值线
        res@cnLevels            =(/20,20.5,21,21.5,22,22.5,23,23.5,24,24.5,25,25.5,26,26.5,27,27.5,28,28.5,29,29.5,30,30.5/)             ;22个值
        res@cnFillColors        =(/24,32,37,44,51,61,70,74,81,97,113,138,145,153,161,166,174,181,186,190,196,199,215/)                 ;23个值,要比上一行多一个值
        plots   = new(5,graphic)                                             ;创建要画图的数量
                plots(0)=gsn_csm_contour_map(wks,R5,res)        ;画图1,gsn_csm_contour_map是一个画填色图的函数,res是指以上的画图设置,wks是前面的一个设置,设置了文件位置,R5是wrf模式资料里前面设置的一个变量,表示月平均温度
                plots(1)=gsn_csm_contour_map(wks,ave1,res)      ;画图2,gsn_csm_contour_map是一个画填色图的函数,res是指以上的画图设置,wks是前面的一个设置,设置了文件位置,ave1是ncep资料里前面设置的一个变量,表示月平均温度
                plots(2)=gsn_csm_contour_map(wks,aave2,res)     ;画图3,gsn_csm_contour_map是一个画填色图的函数,res是指以上的画图设置,wks是前面的一个设置,设置了文件位置,aave2是ec资料里前面设置的一个变量,表示月平均温度
                plots(3)=gsn_csm_contour_map(wks,ave3,res)      ;同上
                plots(4)=gsn_csm_contour_map(wks,R,res)         ;同上
;---Arrays to hold text annotation ids(这里开始是为了要设置每张图左上角的字母)
        txid_tl                                    = new(5,graphic)        
        amid_tl                                  = new(5,graphic)
        txres                                     = True
        txres@txPerimOn                   = False                  
        txres@txFontHeightF              = 0.04          ;图内字体的大小
        txres@txFontColor                  ="white"        ;图内字体颜色
        txres@txFontThicknessF          =37             ;字的粗细,可以不设置
;---Top left string
        amres_tl                               = True           
        amres_tl@amParallelPosF       = -0.49         ;调整字符距离左边的远近
        amres_tl@amOrthogonalPosF = -0.49         ;调整字符距离离上边的远近.
        amres_tl@amJust                  = "TopLeft"     ;字符的位置,这里设置的是左上角
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        do i=0,4                                        ;写循环是为了让字符可以按顺序写a,b,c,d,e,这里设置0,4是因为这里有五张图,需要进行5次循环
;---Create text strings
                        if (i .eq. 0) then
                                tl_label = "(a)"
                        else if (i .eq. 1) then
                                tl_label = "(b)"
                        else if (i .eq. 2) then
                                tl_label = "(c)"
                        else if (i .eq. 3) then
                                tl_label = "(d)"
                        else if (i .eq. 4) then
                                tl_label = "(e)"
                        end if
                        end if
                        end if
                        end if
                        end if
                        txid_tl(i) = gsn_create_text(wks, tl_label, txres)                     ;将字符写进图里
                        amid_tl(i) = gsn_add_annotation(plots(i), txid_tl(i), amres_tl)  
                end do

        pres                               = True
        pres@gsnMaximize          = True
        pres@gsnPanelLabelBar   = True                 ;是否显示总的色板
        pres@pmLabelBarWidthF   = 0.8                 ;色板宽度
        pres@lbLabelFontHeightF = 0.008               ;坐标中数字的大小
        gsn_panel(wks,plots,(/3,2/),pres)               ;这里是图摆放的方式。(3,2)是指五张图排成3行2列
end






密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2019-6-5 22:28:32 | 显示全部楼层
好厉害,在研究怎么能把这个帖子收藏了呢
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-6-6 15:02:49 | 显示全部楼层
感谢楼主 非常优秀 很有用了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-8-3 12:01:26 | 显示全部楼层
非常棒 感谢楼主
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-1-24 16:25:29 | 显示全部楼层
谢谢楼主
一定多多学习
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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