ncl空间分布图代码

veroniaca

如果想要画出类似下面的图片，您可以好好看看如下代码：（黄色背景都是数据的处理，可以直接不看。绿色字是无关紧要的口令，可直接省略。其他均是绘图的设置）这个代码出来的图片，就是我下面贴出来的这张。
·
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

CDH

好厉害，在研究怎么能把这个帖子收藏了呢

supersupergirl

感谢楼主 非常优秀 很有用了

supersupergirl

非常棒 感谢楼主

Growing

谢谢楼主
一定多多学习