NCL站点资料画中国区域气温散点图及分布图
昨天发了按书上画的,我感觉我学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。
[email protected]_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"
[email protected]_name = "lon"
[email protected] = "degrees-east"
lon_av&lon = lon_av
lat_av!0 = "lat"
[email protected]_name = "lat"
[email protected] = "degrees_north"
lat_av&lat = lat_av
olon!0 = "lon"
[email protected]_name = "lon"
[email protected] = "degrees-east"
olon&lon = olon
olat!0 = "lat"
[email protected]_name = "lat"
[email protected] = "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
[email protected]=False ;如果设置为真,则循环点被加入数据,如果数据不是循环的,就设置为假就可以。
[email protected] = "Ncarg4_1";网上的那个代码里没有这句,害我折腾了好久才明白
[email protected]="Earth..4" ;中国地图包含在这个叫Earth..4的地图库里
[email protected]=True
[email protected]=(/"China:states","*"/)
[email protected]=2.0 ;这两行是为了加粗边界和国界线
[email protected]=2.0
[email protected]=15.0
[email protected]=55.0
[email protected]=72
[email protected]=135.0
[email protected] = "Ncarg4_1"
[email protected] = True ;使能填充覆盖
[email protected] = (/"China:states","*"/)
[email protected]=0
[email protected]=0
[email protected]=True;画填充图
[email protected]=False;不画等值线
[email protected]=False;不要等值线上的标签
[email protected]="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。
[email protected]_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
[email protected]=False
[email protected]="Always"
[email protected]=15.0
[email protected]=55.0
[email protected]=72
[email protected]=135.0
[email protected]="1951 January China average temperature"
[email protected] = "Ncarg4_1";这一步和下一步必须要,否则加载中国地图的时候会出错(找不到地图库)
[email protected]="Earth..4";这步加上步再加下面那个China:state和*就可以画出中国轮廓边界了
[email protected]=True
[email protected]=(/"China:states","*"/);在这个地图库里我们绘制中国和*的区域边界,China:state里居然没有*!不能忍,要包括进来!
[email protected]=2.0 ;这两行是为了加粗边界和国界线
[email protected]=2.0
;我看网上那个最多人转载的那个站点画气温分布等值线图,它用的兰伯特投影,我画了一下,左右纬度不对称,很奇怪,就放弃使用了。
map=gsn_csm_map(wks,mpres)
gsres=True
[email protected]=16
do i=0,num_markers-1
if (.not.ismissing(lat_new(i,0))) then
[email protected]=colors(i)
[email protected]=i+1
gsn_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres)
end if
end do
frame(wks);
end
画的是散点分布图,图如下:
初步入门了NCL,真是令人高兴!还要继续努力!