NCL 用站点数据作图，有3个waring，而且无等值面，提示无有效数据

gycsjz

gycsjz

脚本警告提示如下：
warning:Argument 0 of the current function or procedure was coerced to the appropriate type and thus will not change if the function or procedure modifies its value
warning:undef: attempting to undefine a key word
warning:Argument 0 of the current function or procedure was coerced to the appropriate type and thus will not change if the function or procedure modifies its value
warning:ContourPlotInitialize: no valid values in scalar field; ContourPlot not possible:[errno=1101]

我的脚本：
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"

load "$GEODIAG_ROOT/geodiag.ncl"

begin
    undef(imagename)         ;接收ncl command传入文件名
    undef(fname)            ;接收数据文件
    undef(coordinateargu);接收ncl comman 传入经纬度范围数组
    undef(flag)             ;
    undef(isnanhai)         ;传入是否显示南海标识 False -否 True - 是
    undef(arraylevels)      ;接收传入色标等级数组
    undef(colormap)             ;接收传入色标数组

    minlat = coordinateargu(0)          ;取得维度最小值

    maxlat = coordinateargu(1)          ;取得维度最大值

    minlon = coordinateargu(2)          ;取得经度最小值

    maxlon = coordinateargu(3)          ;取得经度最大值

    lines = asciiread(fname,-1,"string")
    lon = stringtofloat(str_get_field(lines(0:),3," "));经度数据
    lat = stringtofloat(str_get_field(lines(0:),4," "));维度数据
    pwv = stringtofloat(str_get_field(lines(0:),5," "));pwv数据

    if (flag .eq. "1") then
    pwv = stringtofloat(str_get_field(lines(0:),6," "))
    else if (flag .eq. "3") then
    pwv = stringtofloat(str_get_field(lines(0:),7," "))
    end if
    end if


    olon = new(66,"float");
    olat = new(40,"float");
    data1 = new((/40,66/),"float")

    do i=0,65
       olon(i) =72+i
    end do
    do l=0,39
       olat(l) = 17+l
    end do


    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

    pwv@_FillValue = 999999.000000  
    rscan = (/10,5,3/)   ;连续的有效半径大小，最大为10，依次递减
    data1 = obj_anal_ic_deprecated(lon,lat,pwv,olon,olat,rscan, False)  ;Creanm插值
; print(data1)
    cmap = RGBtoCmap(colormap)                      ; 取得自定义的colormap颜色
    wtype = "png"
    wtype@wkWidth = 600
    wtype@wkHeight = 419
    wks = gsn_open_wks(wtype, imagename)            ;创建工作空间
    gsn_define_colormap(wks,cmap)

;   设置地图
;>--------------------------------------------<
    res=True
    res@gsnMaximize           = True
    res@gsnPaperOrientation   = "portrait"
    res@isUseLowResolution = False
    res@isShowProvince = True
    res@isShowCity = False
    res@cityThickness = 2.0
    if (isnanhai .eq. "False") then
    res@isShowSouthChinaSea = False
    else
    res@isShowSouthChinaSea = True
    end if
   
    res@isAddMask = True
    res@isShowRivers = True
    res@riverColor = "black"
    res@riverThickness = 1.0
    res@boundaryThickness = 2.0
    res@southChinaSeaMapWidth = 0.1

;<---------设置投影------------------->
    res@gsnAddCyclic  = False
    res@mpProjection = "Mercator";地图投影方式,墨卡托
    res@mpLambertMeridianF = 110.0
    res@mpLimitMode = "LatLon"
    res@mpLambertParallel1F = .001      ;Default: .001
    res@mpLambertParallel2F = 89.999    ;Default: 89.999
    res@mpOceanFillColor = "white"     ;用白色填充海洋
    res@mpInlandWaterFillColor = "white"  ;用白色填充内陆湖水
    res@mpLandFillColor       = "white" ;用白色填充陆地
    res@mpInlandWaterFillColor = "white"

    res@mpMinLatF             =  minlat
    res@mpMaxLatF             =  maxlat
    res@mpMinLonF             =  minlon
    res@mpMaxLonF             =  maxlon

;<----------设置网格-------------------->
;   res@mpGridAndLimbOn        = True
;   res@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
;   res@mpGridLineDashPattern  = 5                  ; lat/lon lines dashed
;   res@mpGridLatSpacingF      = 5
;   res@mpGridLonSpacingF      = 5
    res@mpGridAndLimbDrawOrder = "PreDraw"          ; Draw grid first
   
    setup_china_map(res) ; <--- set the resource for China map


    ; --------------------------------------------------------------------------
    ; draw any data you want here

    res@cnFillDrawOrder = "PreDraw"
    res@cnFillOn        = True  ;等值线填充
    res@cnLinesOn       = False ;等值线不显示
    res@cnLineLabelsOn  = False ;等值线标签不显示
;   res@cnLevelSpacingF = 5.   ;等值线间距
    res@cnSpanFillPalette = True

    res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
    res@cnLevels             = arraylevels

    res@lbLabelAutoStride = True
    res@lbLabelBarOn    = True
    res@lbAutoManage = False
    res@pmLabelBarDisplayMode = "Always"
    res@lbBoxLinesOn = False
    res@pmLabelBarHeightF = 0.06
;    res@pmLabelBarWidthF = 0.4
    res@lbPerimOn = False

    res@gsnLeftString = ""  ; 关闭左上角注释信息
    res@gsnRightString = ""  ;关闭右上角注释信息

    res@pmTickMarkDisplayMode = "Always" ;刻度线
    res@tmXTOn = False                   ;关闭X顶部刻度线
    res@tmYROn = False                   ;关闭Y右部刻度线
;    res@tmXBLabelFontHeightF = 0.06
    res@tmXBMajorLengthF = 0.005
;   res@gsnMinorLonSpacing = 5


;   对图片的设置
     
    res@vpXF = 0.8    ;左边距
    res@vpYF = 0.95 ;上边距
    res@vpWidthF = 1.0
    res@vpHeightF = 0.8
    plot = gsn_csm_contour_map(wks, data1, res)
   
    ; --------------------------------------------------------------------------
    ; attach China map on the base plot
    attach_china_map(wks, plot)

    draw(plot)
    frame(wks)
end
exit

所用数据格式为：

12345 hhhh 36.905551304 100.99010042 10.6 99999 99999
12345 aaaa 27.693595905 117.41039069 41.5 99999 99999


琢磨不通了。测试发现插值后data1的值都为@FillValue设置的999999。不知道为什么。忘大侠们指点

言深深

帮楼主合并了此2贴

longlivehj

    pwv = stringtofloat(str_get_field(lines(0:),5," "));pwv数据

    if (flag .eq. "1") then
    pwv = stringtofloat(str_get_field(lines(0:),6," "))
    else if (flag .eq. "3") then
    pwv = stringtofloat(str_get_field(lines(0:),7," "))
    end if
    end if

感觉你上面这段逻辑有问题。想问一下，flag的意思？另外，第6、7列是不是都是99999啊？

gycsjz

longlivehj 发表于 2014-8-25 18:09
pwv = stringtofloat(str_get_field(lines(0:),5," "));pwv数据

    if (flag .eq. "1") then

flag是传入的字符串0.1.3根据不同情况取得不同列的数据。我传入0取得就是第5列。1就是第6列。3就是第7列这个意思。6.7列是5个9嘛。不是缺省值6个9

longlivehj

gycsjz 发表于 2014-8-25 16:37
脚本警告提示如下：
warning:Argument 0 of the current function or procedure was coerced to the appro ...

12345 hhhh 36.905551304 100.99010042 10.6 99999 99999
12345 aaaa 27.693595905 117.41039069 41.5 99999 99999

从测试数据上看，第3列是纬度，第4列是经度。

lon = stringtofloat(str_get_field(lines(0:),3," "));经度数据
lat = stringtofloat(str_get_field(lines(0:),4," "));维度数据

gycsjz

longlivehj 发表于 2014-8-26 09:05
12345 hhhh 36.905551304 100.99010042 10.6 99999 99999
12345 aaaa 27.693595905 117.41039069 41.5  ...

我了个去。还看半天写的程序。原来库中的经纬度入反了。谢过谢过

gycsjz

longlivehj 发表于 2014-8-26 09:05
12345 hhhh 36.905551304 100.99010042 10.6 99999 99999
12345 aaaa 27.693595905 117.41039069 41.5  ...

谢了。确实是这问题

longlivehj

gycsjz 发表于 2014-8-26 09:48
谢了。确实是这问题

呵呵，不客气！

另外，
    do i=0,65
       olon(i) =72+i
    end do
可以写成
    olon = 72 + ispan(0, 65, 1)

最好用obj_anal_ic_Wrap替代 obj_anal_ic_deprecated，后者实现不完整，注意最后的一个单词deprecated！

gycsjz

longlivehj 发表于 2014-8-26 09:53
呵呵，不客气！

另外，

多谢指导。有时间研磨研磨