NCL 画台风路径图问题

指甲钳 · 发表于 2013-11-13 14:47:11

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由指甲钳于 2013-11-13 14:48 编辑

想用NCL画台风路径图，看到网上一篇很详细的步骤，http://blog.sina.com.cn/s/blog_a134572301014xft.html，用的是JTWC的数据，TXT文本格式的，然后我想用自己的数据画。但一开始就出错了。但想不通为什么出错啊，想请教一下大家。

我的数据内容：第一列时间，第二列经度，第三列维度，第四列气压，第五列风速

2013110408, 152.2, 6.2, 1000,       18,
2013110414, 150.2, 6.3, 995,       20,
2013110420, 149, 6.4, 995,       20,
2013110502, 147.6, 6.5, 985,       25,
2013110508, 146.1, 6.6, 980,       30,
2013110514, 144.6, 6.6, 975,       33,
2013110520, 143.1, 6.7, 970,       35,
2013110602, 141.3, 7.2, 955,       42,
2013110608, 139.7, 7.4, 930,       55,
2013110614, 137.9, 7.6, 925,       58,
2013110620, 136.2, 7.9, 915,       62,
2013110702, 134.3, 8.2, 905,       65,
2013110708, 132.8, 8.7, 900,       68,
2013110714, 131.1, 9.3, 895,       70,
2013110720, 129.1, 10.2, 890,       75,
2013110802, 127, 10.6, 890,       75,
2013110808, 124.8, 11, 892,       72,
2013110811, 123.7, 11.3, 920,       60,
2013110814, 122.6, 11.4, 940,       50,
2013110817, 121.6, 11.5, 940,       50,
2013110820, 120.7, 11.8, 940,       50,
2013110823, 119.2, 12.2, 940,       50,
2013110902, 118.2, 12.5, 945,       48,
2013110905, 117.4, 12.5, 950,       45,
2013110908, 116.8, 12.5, 950,       45,
2013110909, 116.5, 12.5, 950,       45,
2013110910, 116.1, 12.6, 950,       45,
2013110911, 115.8, 12.8, 950,       45,
2013110912, 115.6, 13, 950, 45       ,
2013110913, 115.3, 13.3, 950, 45,
2013110914, 115, 13.5, 950, 45       ,
2013110915, 114.7, 13.6, 950, 45,
2013110916, 114.5, 13.7, 950, 45,
2013110917, 114.3, 13.8, 950, 45,
2013110918, 114.1, 14, 950, 45       ,
2013110919, 113.8, 14.2, 950, 45,
2013110920, 113.5, 14.5, 950, 45,
2013110921, 113.1, 14.6, 950, 45,
2013110922, 112.6, 14.8, 950, 45,
2013110923, 112.3, 14.9, 955, 42,
2013111000, 112, 15.2, 955, 42       ,
2013111001, 111.8, 15.4, 955, 42,
2013111002, 111.6, 15.4, 955, 42,
2013111003, 111.4, 15.6, 955, 42,
2013111004, 111.1, 15.8, 955, 42,
2013111005, 110.8, 15.9, 955, 42,
2013111006, 110.7, 16.1, 955, 42,
2013111007, 110.4, 16.2, 955, 42,
2013111008, 110.2, 16.4, 955, 42,
2013111009, 110, 16.7, 955, 42,
2013111010, 109.8, 16.9, 955, 42,
2013111011, 109.5, 17.1, 955, 42,
2013111012, 109.4, 17.3, 955, 42,
2013111013, 109.2, 17.5, 955, 42,
2013111014, 109,17.8, 955, 42,
2013111015, 108.8, 18.1, 955, 42,
2013111016, 108.5, 18.4, 955, 42,
2013111017, 108.3, 18.6, 955, 42,
2013111018, 108.1, 18.7, 955, 42,
2013111019, 108, 18.9, 955, 42       ,
2013111020, 107.8, 19, 955, 42       ,
2013111021, 107.7, 19.3, 960, 40,
2013111022, 107.6, 19.6, 960, 40,
2013111023, 107.6, 19.8, 965, 38,
2013111100, 107.6, 20, 965, 38,
2013111101, 107.5, 20.3, 965, 38,
2013111102, 107.4, 20.5, 965, 38,
2013111103, 107.4, 20.7, 965, 38,
2013111104, 107.3, 20.9, 965, 38,
2013111105, 107.3, 21.1, 965, 38,
2013111106, 107.3, 21.3, 970, 35,
2013111107, 107.3, 21.5, 970, 35,
2013111108, 107.3, 21.6, 975, 33,
2013111109, 107.3, 21.8, 975, 33,
2013111110, 107.4, 21.9, 980, 30,
2013111111, 107.6, 22, 985, 28,
2013111112, 107.8, 22, 985, 28,
2013111113, 107.9, 22.2, 985, 28,
2013111114, 108, 22.4, 988, 25,
2013111115, 108.1, 22.4, 988, 25,
2013111116, 108.2, 22.4, 990, 23,
2013111117, 108.3, 22.5, 992, 20,
2013111118, 108.4, 22.5, 992, 20,
2013111119, 108.5, 22.5, 994, 18,
2013111120, 108.6, 22.5, 996, 16,


脚本：
;********************************************************
; Load NCL scripts
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

  fiTY = "e:/sql_data/haiyan.txt"

; 获取文本文件的行数，相应的还有numAsciiCol函数用于获取列数
  nrow = numAsciiRow(fiTY)

  YYYYMMDDHH = new(nrow, "string")
  lat = new(nrow, "float")
  lon = new(nrow, "float")
  vmax = new(nrow, "integer")
  mslp = new(nrow, "integer")

  data = asciiread(fiTY, -1, "string")
  YYYYMMDDHH = str_get_field(data, 1, ",")
  lat = stringtofloat(str_get_field(data,3, ","))
  lon = stringtofloat(str_get_field(data, 1, ","))
  vmax = stringtoint(str_get_field(data, 5, ","))
  mslp = stringtoint(str_get_field(data, 4, ","))

  DateChar = stringtochar(YYYYMMDDHH)
  MM = chartostring(DateChar(:,5:6))
  DD = chartostring(DateChar(:,7:8))
  HH = chartostring(DateChar(:,9:10))

  HHi = mod(stringtoint(YYYYMMDDHH), 100)
  DDi = mod(stringtoint(YYYYMMDDHH)/100, 100)
  MMi = mod(stringtoint(YYYYMMDDHH)/10000, 100)

; plot

  wks = gsn_open_wks("png", "haiyan")
  gsn_define_colormap(wks,"rainbow")

  res = True

  res@gsnDraw = False
  res@gsnFrame = False

; 设置底图边界
  res@mpMinLatF = 10
  res@mpMaxLatF = 50
  res@mpMinLonF = 110
  res@mpMaxLonF = 160

; 设置海陆填充颜色，这里的颜色选取rainbow色表中的第160号颜色
; 海洋默认是白色
  res@mpLandFillColor = 160
; res@mpOceanFillColor = "white"

; 绘制洋面经纬网格
  res@mpGridAndLimbOn = "True"
  res@mpGridMaskMode = "MaskNotOcean"
  res@mpGridLineDashPattern = 15
  res@mpGridSpacingF = 5.0

; 绘制国界
  res@mpOutlineOn = True
  res@mpOutlineBoundarySets = "National"

; 绘制省界
  res@mpDataBaseVersion = "MediumRes"
  res@mpDataSetName = "Earth..4"
  res@mpOutlineSpecifiers = "China:States"

  plot = gsn_csm_map_ce(wks,res)

;================
; 按照vmax(单位：节 knot，海里/小时)变量提供的风速大小对台风强度进行区分，
; 并在绘制时用不同颜色标出
; 0~33.5 热带低压
; 34~63 热带风暴
; 64~80 台风
; 81~99 强台风
; >100 超强台风

  colours = (/3, 20, 60, 120, 180/)

  resDot = True
  resLine = True

  dumDot= new(nrow, graphic)
  dumLine = new(nrow, graphic)

; 绘制线
  resLine@gsLineThicknessF = 3

  do i = 0, nrow-2
xx = (/ lon(i), lon(i+1)/)
yy = (/ lat(i), lat(i+1)/)

if (vmax(i).le.33) then
   resLine@gsLineColor = colours(0)
end if
if (vmax(i) .ge. 34 .and. vmax(i).le.63) then
   resLine@gsLineColor = colours(1)
end if
if (vmax(i).ge.64 .and. vmax(i) .le. 80) then
   resLine@gsLineColor = colours(2)
end if
if (vmax(i).ge.81 .and. vmax(i) .lt. 99) then
   resLine@gsLineColor = colours(3)
end if
if (vmax(i).ge.100) then
   resLine@gsLineColor = colours(4)
end if

dumLine(i) = gsn_add_polyline(wks, plot, xx, yy, resLine)
  end do

; 绘制00时的点
  resDot@gsMarkerIndex = 1
  resDot@gsMarkerSizeF = 0.02
  resDot@gsMarkerColor = "black"

  do i = 0, nrow-1

if (HH(i) .eq. "00") then
   dumDot(i) = gsn_add_polymarker(wks, plot, lon(i), lat(i), resDot)
end if

  end do

; 绘制图例
  resLg = True

  resLg@lgItemType = "MarkLines"

  resLg@lgMonoMarkerIndex = True
  resLg@lgMarkerColors = colours
  resLg@lgMarkerIndex = 1
  resLg@lgMarkerSizeF = 0.02

  resLg@lgMonoDashIndex = True
  resLg@lgDashIndex = 0
  resLg@lgLineColors = colours
  resLg@lgLineThicknessF = 3

  resLg@vpWidthF = 0.14
  resLg@vpHeightF = 0.13

  resLg@lgPerimFill = 0
  resLg@lgPerimFillColor = "Background"

  resLg@lgLabelFontHeightF = 0.08

  resLg@lgTitleFontHeightF = 0.015
  resLg@lgTitleString = "Best Track"

  lbid = gsn_create_legend(wks, 5, (/" 33kt or less",  " 34 to 63kt", " 64 to 80kt", " 81 to 99kt",  " 100 or more"/), resLg)

; 将图例放置在图中
  amres = True
  amres@amParallelPosF = -0.3
  amres@amOrthogonalPosF = -0.3
  dumLg = gsn_add_annotation(plot, lbid, amres)

; 标注00时的日期
  dumDate = new(nrow,graphic)

  resTx = True
  resTx@txFontHeightF = 0.012
  resTx@txFontColor = "black"
  resTx@txJust = "BottomLeft"

  do i = 1, nrow-1
if (HH(i) .ne. "00" ) then
   continue
end if
dumDate(i) = gsn_add_text(wks,plot, MM(i)+DD(i), lon(i)+0.5, lat(i), resTx)
  end do

  draw(wks)
  frame(wks)

end

基本上没有改动，但为什么出现提示错误：
  YYYYMMDDHH = str_get_field(data, 1, ",")这个语句，左右的维数不匹配？
Dimension sizes of left hand side and right hand side of assignment do not match

二爷的屁股 · 发表于 2013-11-13 16:01:04

帮顶啊啊啊啊啊啊啊啊

指甲钳 · 发表于 2013-11-13 16:18:46

二爷的屁股发表于 2013-11-13 16:01

登录/注册后可看大图

帮顶啊啊啊啊啊啊啊啊

谢谢啊啊啊啊啊啊啊

尽头的尽头 · 发表于 2013-11-14 14:16:19

你之前这句不是YYYYMMDDHH = new(nrow, "string")已经定义YYYYMMDDHH的维数跟你的数据一样
YYYYMMDDHH = str_get_field(data, 1, ",")这句获取数据是一维的，当然不匹配了

指甲钳 · 发表于 2013-11-14 14:41:05

尽头的尽头发表于 2013-11-14 14:16

登录/注册后可看大图

你之前这句不是YYYYMMDDHH = new(nrow, "string")已经定义YYYYMMDDHH的维数跟你的数据一样
YYYYMMDDHH = s ...

这样啊。谢谢啊！

尽头的尽头 · 发表于 2013-11-14 14:44:43

指甲钳发表于 2013-11-14 14:41

登录/注册后可看大图

这样啊。谢谢啊！

不客气~~

尽头的尽头 · 发表于 2013-11-14 15:16:03

我吧你的数据复制下来做的时候，也发现这个问题，你的数据明明有85行，但是用ncl的numAsciiRow获取出来只有84，所以不匹配了

指甲钳 · 发表于 2013-11-14 17:16:02

尽头的尽头发表于 2013-11-14 15:16

登录/注册后可看大图

我吧你的数据复制下来做的时候，也发现这个问题，你的数据明明有85行，但是用ncl的numAsciiRow获取出来只有 ...

我下载过原文的数据，数据的变量比我多，但格式应该是一样的，但原文的数据就可以画出图来，会不会是以为他的列数从0开始，不是从1开始？但这样的话原文应该画不出来啊，我看过也有网友出过一样的问题，但未找找到答案，所以上来问问。

尽头的尽头 · 发表于 2013-11-14 17:53:28

指甲钳发表于 2013-11-14 17:16

登录/注册后可看大图

我下载过原文的数据，数据的变量比我多，但格式应该是一样的，但原文的数据就可以画出图来，会不会是以为 ...

从什么开始无所谓，关键是少了一行，这就是问题，不用他那种方法倒是可以画图

指甲钳 · 发表于 2013-11-14 18:59:05

尽头的尽头发表于 2013-11-14 17:53

登录/注册后可看大图

从什么开始无所谓，关键是少了一行，这就是问题，不用他那种方法倒是可以画图

怎么画啊？求指教！！

		自动登录	找回密码
密码			立即注册

[作图] NCL 画台风路径图问题

登录后查看更多精彩内容~