参考：http://bbs.06climate.com/forum.php?mod=viewthread&tid=11417&extra=&page=1

1、读取站点数据

lines = asciiread("e:\a.txt",-1,"string") ;站点文件绝对路径e:\a.txt，为ascii文件

;即打开后可显示正常文字的txt、dat等文件

;-1表示不考虑维数

;string表示读入的数据视作字符串类型

;得到的lines为n*1维的数组，n为行数

lon = stringtofloat(str_get_field(lines(0),2," ")) ;读取第一行第2个数据，并转换为浮点型

2、站点插值到网格点

a) 定义经纬度变量

olon = new(51,"float") ;定义横坐标（经度）有51个格点，值均为浮点型

olat = new(31,"float") ;定义纵坐标（纬度）有31个格点

 

 

do i=0.50 ;对新建的经纬度赋值

olon(j) = 110+j*0.25 ;经度范围为110°~122.5°E，格点精度为0.25°

end do

 

do i=0,30

olat(j) = 10+j*0.25 ;纬度范围为10°~17.5°N，格点精度为0.25°

end do

 

olon!0          = "lon" ;经度为一维变量，第一维名称为lon

olon@long_name  = "lon" ;告诉ncl第一维代表经度

olon@units      = "degrees-east" ;以向东为正方向

olon&lon        = olon ;经度的值取olon

 

olat!0          = "lat"

olat@long_name  = "lat"

olat@units      = "degrees_north"

olat&lat        = olat

b) 定义与站点对应的格点变量

grid_T = new((/31,51/),"float") ;定义格点场为31*51的二维变量

 

c) 插值

T@_FillValue = 9999.0 ;定义原变量T中的缺测值

;凡是记录值等于9999.0的点均忽略

rscan = (/.5 , .3 , .15 , .05/) ;定义影响半径数组

;数值单位为度，数组内数值依次递减
;数值越大，插值得到的结果越平滑

grid_T = obj_anal_ic_deprecated(lon,lat,T,olon,olat,rscan,False)

;利用Cressman插值，得到插值后的格点场grid_T

 

3、设置resource，绘制grid_T

另：标记站点位置的画法

1. 用map = gsn_csm_map(wks,res)画出地图

2. 读入站点文件的经纬度lon0,lat0

3. 设置站点的resource，令其为mkres

mkres = True

mkres@gsMarkerIndex = 16 ;标记的种类，16号为填色圆圈

mkres@gsMarkerSizeF = 0.003 ;设定标记的大小

mkres@gsMarkerColor = 230 ;设定标记颜色为colormap的第231个颜色

 

gsn_polymarker(wks,map,lon0,lat0,mkres) ;绘制

http://bbs.06climate.com/data/attachment/album/201403/06/183018lsdtt8yonxnytlld.png

Marker Table