昨天发了按书上画的，我感觉我学NCL最大问题在于，毫无基础的情况下，想像学MATLAB那样，学了基础语法以后，直接看网上具有某特定功能的程序，然后遇到函数再查相关问题，然而问题在于MATLAB基数大，NCL用的人少，遇到的问题少，因此不容易查到和你遇到相同error或者说相关提问很少，对debug困难很大。不过NCL官网上各种函数、各种错误介绍的都很全就是啦。我之前借的学姐的书，我买的书今天才到（《NCL图形分析语言入门到精通》）。我感觉学编程语言应该先从书上学习最简单的例子，一开始就仿制那些复杂的、好看的操作不利于了解函数、句柄。我在气象家园网站上仿作别人的，结果太复杂了，报的错我没法debug，就束手无策。（实际证明原代码少一行！）而按书上做的就非常清楚。

今天我把我读的某一个txt文件。绘制单要素的地图散点图的代码贴上来，实际上也是根据昨天最原始的那个编制的，应该不困难。其中我参考了网上代码：http://bbs.06climate.com/forum.php?mod=viewthread&tid=11417 最NCL做站点资料不妨看看这个。代码如下：

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"
begin
    
    file_path="/mnt/c/users/hong/desktop/1.txt"
    data=asciiread(file_path, (/4216,13/), "float")
    station=data(:,0)       ;读入站点号
    lat=data(:,1)           ;读入纬度
    lon=data(:,2)           ;读入经度
    day=data(:,6)           ;由于我选的就是1951年1月的，所以只有日期存在区别。
    temp=data(:,7)          ;NCL是从0开始计数，因此第8列索引是7。
    temp@_FillValue=32766   ;在数据说明里说了气温的缺测值是32766
     ;由于在资料里是整数性，需要对其进行一下转换，其中经纬度要从60进制转为100进制
    lat_a=round(lat*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100
    lon_a=round(lon*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100
    temp_a=temp*0.1
    temp_av=new(136,"float")
    lon_av=new(136,"float")
    lat_av=new(136,"float")
    do i=0,135
        temp_av(i)=dim_avg_n(temp_a(31*i:31*i+30),0)
        lon_av(i)=lon_a(31*i)
        lat_av(i)=lat_a(31*i)
    end do
    temp_max=max(temp_av)
    temp_min=min(temp_av)
    olon=new(64,"float")    ;中国经度范围73°33′E-135°05′E，这里我设置经度70°-140°
    olat=new(41,"float")    ;中国纬度范围3°51′N-53°33′N，这里我设置纬度0°-55°
    data1=new((/41,64/),"float")
    do i=0,63
        olon(i)=72.0+i
    end do
    do l=0,40
        olat(l)=l+15
    end do
    ;要设置经纬度的单位和名称，否则会报错
     lon_av!0          = "lon"
     lon_av@long_name  = "lon"
     lon_av@units      = "degrees-east"
     lon_av&lon        = lon_av
     lat_av!0          = "lat"
     lat_av@long_name  = "lat"
     lat_av@units      = "degrees_north"
     lat_av&lat        = lat_av
     olon!0          = "lon"
     olon@long_name  = "lon"
     olon@units      = "degrees-east"
     olon&lon        = olon
     olat!0          = "lat"
     olat@long_name  = "lat"
     olat@units      = "degrees_north"
     olat&lat        = olat
     ;进行插值操作
    rscan=(/20,10,5/);网上是10 5 3 我这里用这个是为了让每个位置有数字，因为我没设置缺测值
    R=obj_anal_ic_deprecated(lon_av,lat_av,temp_av,olon,olat,rscan, False)
    ;下面在开工作空间之前，设置了colormap，我直接复制过来了。
    cmap = (/                     \
            (/ 255./255, 255./255, 255./255 /),    \  ; 0 - White background.
            (/ 0./255  , 0./255  , 0./255   /),         \  ; 1 - Black foreground.
            (/ 255./255, 0./255  , 0./255   /),        \  ; 2 - Red.
            (/ 0./255  , 0./255  , 255./255 /),        \  ; 3 - Blue.
            (/ 164./255, 244./255, 131./255 /),      \  ; 4 - Ocean Blue.
            (/ 0./255  , 0./255  , 255./255   /),         \  ; 5 - Bar 1
         (/ 0./255  , 153./255, 255./255 /),      \  ; 6 - Bar 2
            (/ 0./255, 153./255, 153./255 /),        \  ; 7 - Bar 3 
            (/ 0./255  , 255./255, 0./255  /),       \  ; 8 - Bar 4   
            (/ 255./255, 255./255  , 102./255 /),        \  ; 9 - Bar 5 
            (/ 255./255, 153./255  , 102./255   /),        \  ; 10 - Bar 6
            (/ 255./255, 0./255  , 255./255   /)        \  ; 11 - Bar 7   
         /)

    wks=gsn_open_wks("png","example_195101")
    gsn_define_colormap(wks,cmap)      ; 
    res=True
    res@gsnAddCyclic=False             ;如果设置为真，则循环点被加入数据，如果数据不是循环的，就设置为假就可以。
    res@mpDataBaseVersion = "Ncarg4_1";网上的那个代码里没有这句，害我折腾了好久才明白
    res@mpDataSetName="Earth..4"       ;中国地图包含在这个叫Earth..4的地图库里
    res@mpOutlineOn=True
    res@mpOutlineSpecifiers=(/"China:states","Taiwan"/)
    res@mpGeophysicalLineThicknessF=2.0 ;这两行是为了加粗边界和国界线
    res@mpNationalLineThicknessF=2.0
    res@mpMinLatF=15.0
    res@mpMaxLatF=55.0
    res@mpMinLonF=72
    res@mpMaxLonF=135.0

    res@mpDataBaseVersion = "Ncarg4_1"
    res@mpAreaMaskingOn = True   ;使能填充覆盖
    res@mpMaskAreaSpecifiers = (/"China:states","Taiwan"/)
    res@mpOceanFillColor=0
    res@mpInlandWaterFillColor=0

    res@cnFillOn=True;画填充图
    res@cnLinesOn=False;不画等值线
    res@cnLineLabelsOn=False;不要等值线上的标签
    res@cnFillDrawOrder="PreDraw";先画填充
    map=gsn_csm_contour_map(wks,R,res)

end

画出的是等值线填充图，图如下：

然后是散点分布图代码：

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"
begin
    
    file_path="/mnt/c/users/hong/desktop/1.txt"
    data=asciiread(file_path, (/4216,13/), "float")
    station=data(:,0)       ;读入站点号
    lat=data(:,1)           ;读入纬度
    lon=data(:,2)           ;读入经度
    day=data(:,6)           ;由于我选的就是1951年1月的，所以只有日期存在区别。
    temp=data(:,7)          ;NCL是从0开始计数，因此第8列索引是7。
    temp@_FillValue=32766   ;在数据说明里说了气温的缺测值是32766
     ;由于在资料里是整数性，需要对其进行一下转换，其中经纬度要从60进制转为100进制
    lat_a=round(lat*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100
    lon_a=round(lon*0.01,0)+(lat*0.01-round(lat*0.01,3))/60*100
    temp_a=temp*0.1
    temp_av=new(136,"float")
    lon_av=new(136,"float")
    lat_av=new(136,"float")
    do i=0,135
        temp_av(i)=dim_avg_n(temp_a(31*i:31*i+30),0)
        lon_av(i)=lon_a(31*i)
        lat_av(i)=lat_a(31*i)
    end do
    R=temp_av
    temp_max=max(temp_av)
    temp_min=min(temp_av)
    itv=(temp_max-temp_min)/5
    arr=(/temp_min+itv,temp_min+2*itv,temp_min+3*itv,temp_min+4*itv/);我把所有温度均匀地用4个节点分配为5份
    colors = (/10,38,56,66,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_av(indexes)
            lon_new(i,0:(npts_range-1))=lon_av(indexes)
        end if
        delete(indexes)
    end do
    wks=gsn_open_wks("png", "scatter_example")
    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="1951 January China average temperature"
    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
    ;我看网上那个最多人转载的那个站点画气温分布等值线图，它用的兰伯特投影，我画了一下，左右纬度不对称，很奇怪，就放弃使用了。
    map=gsn_csm_map(wks,mpres)
    gsres=True
    gsres@gsMarkerIndex=16
    do i=0,num_markers-1
        if (.not.ismissing(lat_new(i,0))) then
            gsres@gsMarkerColor=colors(i)
            gsres@gsMarkerThicknessF=i+1
            gsn_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres)
        end if
    end do
    frame(wks);
end

画的是散点分布图，图如下：

初步入门了NCL，真是令人高兴！还要继续努力！