- 积分
- 1655
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-4-12
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2014-4-26 15:09:39
|
显示全部楼层
longlive大神，我自己编了个画闪电的程序，能帮我看看不。感觉我对字符的理解有误~~
这是我的程序的思路，把闪电发生的经度，维度，强度弄成三列放到txt文件里，一共2163行。然后一行一行读取，用gsn_add_text在地图上特定的位置画“+”“-”，强度是正就画正号，负的就画符号。
提示一下：下面程序中中间部分读取冰粒子重要目的是获取模式的经纬度范围，为了更好的对比。这些时间和层数其实没有冰粒子的，画不出来的，就是为了获取底图
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/wrf/WRFUserARW.ncl"
begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
; We generate plots, but what kind do we prefer?
type = "pdf"
; type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"lightning")
gsn_define_colormap(wks,"precip_11lev")
a = addfile("/home/Huanglei/data/d032"+".nc","r")
res = True
res@MainTitle = "REAL-TIME WRF"
; res@gsnDraw = False
;res@gsnFrame = False
mpres = True ; Map resources
mpres@mpOutlineOn = False ; Turn off map outlines
mpres@mpFillOn = False ; Turn off map fill
mpres@mpGridAndLimbOn = True
pltres = True ; Plot resources
pltres@PanelPlot = True ; Tells wrf_map_overlays not to remove overlays
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
it =1 ; TIME LOOP
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
print(it)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
if(isfilevar(a,"QICE"))
qi = wrf_user_getvar(a,"QICE",it)
qi = qi*1000.
qi@units = "g/kg"
end if
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
level = 1 ; LOOP OVER LEVELS
opts = res
opts@cnFillOn = True
opts@gsnSpreadColors = False
opts@ContourParameters = (/ 0.02, 0.5, 0.05 /)
opts@gsnDraw = False
opts@gsnFrame = False
if (isvar("qi"))
qis = qi(level,:,:)
contour = wrf_contour(a,wks,qis,opts)
plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
delete(contour)
end if
;>============================================================<
; add China map
;>------------------------------------------------------------<
shp_name1 = "/home/Huanglei/map/China/diquJie_polyline.shp"
lnres = True
lnres@gsLineColor = "gray25"
lnres@gsLineThicknessF = 0.5
id = gsn_add_shapefile_polylines(wks,plot,shp_name1,lnres)
shp_name2 = "/home/Huanglei/map/China/cnmap/cnhimap.shp"
prres=True
prres@gsLineThicknessF = 2.0
prres@gsLineColor = "black"
plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)
ascii_filename = "/home/Huanglei/data/7241600_20.txt"
seismic = asciiread(ascii_filename,(/2163,3/),"float")
x = seismic(:,0) ; Column 1 of file contains X values.
y = seismic(:,1) ; Column 2 of file contains Y values.
z = seismic(:,2) ; Column 3 of file contains Z values.
do lt=0, 2162
lon=x(lt)
lat=y(lt)
if z(lt)>0
str="+"
end if
if z(lt)<0
str="-"
end if
txres2 = True
txres2@txFont = 2
txres2@txFontHeightF =0.01
txres2@txFontColor = "Red"
txdum1 =gsn_add_text(wks, plot, str, lon,lat, txres2)
draw(plot) ; This will draw the map and the shapefile outlines.
frame(wks)
delete(opts)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end do
end
提示我的if有错，对于这个我没怎么看懂
fatal:Conditional statements (if and do while) require SCALAR logical values, see all and any functions
fatal:["Execute.c":8567]:Execute: Error occurred at or near line 105 in file /home/Huanglei/shandian.ncl
|
|
窥视卡
雷达卡
发表于 2014-4-23 18:56:38
显身卡
，是的是的！
发表于 2014-4-27 15:19:35