Hi,
I try to draw 5 figures on one page , and meet this error :
fatal:Invalid plot ID=34 passed to NhlGetBB
the data used in this script had been transfered on ftp.
Can anyone help me?
Thanks.
JinQ
------------------------------------------------------------------------
ftp> put wrfinput_d01
ftp> put obs_gts_2010-09-18_12_00_00.3DVAR
------------------------------------------------------------------------
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"
;***********************************************************************
; user-defined settings
;***********************************************************************
date = "2010091812"
obfile = "obs_gts_2010-09-18_12_00_00.3DVAR" ; need this data
fgfile = "./wrfinput_d01"
out_type = "ps"
plotname = "./FIG_1p:obs_loc"+date
proc_synop = True
proc_metar = True
proc_ship = True
proc_buoy = True
proc_satem = True
proc_satob = True
proc_gpspw = True
proc_sound = True
proc_airep = True
proc_pilot = True
proc_profl = True
proc_qscat = True
proc_bogus = True
proc_airsr = True
proc_gpsrf = True
proc_radar = True
;***********************************************************************
; end of user-defined settings
;***********************************************************************
if ( isfilepresent(fgfile) ) then
wrf_file=addfile(fgfile+".nc","r")
if ( wrf_file@MAP_PROJ .eq. 1 ) then
mapproj = "LambertConformal"
truelat1 = wrf_file@TRUELAT1
truelat2 = wrf_file@TRUELAT2
clon = wrf_file@STAND_LON
end if
if ( wrf_file@MAP_PROJ .eq. 2 ) then
mapproj = "Stereographic"
truelat1 = wrf_file@TRUELAT1
truelat2 = wrf_file@TRUELAT2
clon = wrf_file@CEN_LON
clat = wrf_file@CEN_LAT
end if
if ( wrf_file@MAP_PROJ .eq. 3 ) then
mapproj = "Mercator"
end if
dsizes = getfiledimsizes(wrf_file)
nx = dsizes(2)
ny = dsizes(3)
xlat=wrf_file->XLAT
xlon=wrf_file->XLONG
lat_ll = xlat(0,0,0)
lat_ur = xlat(0,ny-1,nx-1)
lon_ll = xlon(0,0,0)
lon_ur = xlon(0,ny-1,nx-1)
else
print("Error: no first guess found for retrieving mapping info")
exit
end if
colormap = "BkBlAqGrYeOrReViWh200"
undef("setmpres")
function setmpres(title:string, str1:string, str2:string)
begin
mpres = True
; mpres@gsnPaperOrientation = "portrait"
; mpres@gsnMaximize = False ; Maximize plot in frame.
mpres@gsnDraw = False
mpres@gsnFrame = False ; Don't advance the frame
mpres@mpProjection = mapproj ; choose projection
if ( mapproj .eq. "LambertConformal" ) then
mpres@mpLambertParallel1F = truelat1 ; two parallels
mpres@mpLambertParallel2F = truelat2
mpres@mpLambertMeridianF = clon ; central meridian
end if
if ( mapproj .eq. "Stereographic" ) then
mpres@mpCenterLatF = clat
mpres@mpCenterLonF = clon
end if
mpres@mpLimitMode = "Corners"
mpres@mpLeftCornerLatF = lat_ll
mpres@mpLeftCornerLonF = lon_ll
mpres@mpRightCornerLatF = lat_ur
mpres@mpRightCornerLonF = lon_ur
mpres@pmTickMarkDisplayMode = "Always"
; mpres@tmYROn = False
; mpres@tmXBOn = False
mpres@tmXTMajorLengthF = 0
mpres@tmYLMajorLengthF = 0
mpres@tmXBMajorLengthF = 0
mpres@tmYRMajorLengthF = 0
mpres@mpDataSetName = "Earth..4"
mpres@mpDataBaseVersion = "MediumRes"
mpres@mpOutlineSpecifiers = (/"China:states","India","Taiwan"/)
mpres@mpFillOn = False
;mpres@mpFillOn = True
;mpres@mpLandFillColor = "gray"
;mpres@tfDoNDCOverlay = True
;mpres@mpUSStateLineColor = "gray"
;mpres@mpNationalLineColor = "gray"
;mpres@mpGeophysicalLineColor = "gray"
mpres@mpGridAndLimbOn = True
mpres@mpGridLineDashPattern = 2
mpres@mpGridLineDashSegLenF = 0.06 ; default 0.15
mpres@mpGridMaskMode="MaskFillArea" ; 擦掉网格线
mpres@mpFillColors = (/"background","White","chocolate1","DeepSkyBlue", "transparent"/)
mpres@mpGeophysicalLineColor = "chartreuse3"
mpres@mpUSStateLineColor = "chartreuse3"
mpres@mpUSStateLineThicknessF=2 ;省线粗细
mpres@mpGridMaskMode="MaskFillArea" ; 擦掉网格线
mpres@mpLimbLineColor = "Red"
mpres@mpNationalLineColor = "Blue"
mpres@mpPerimLineColor = "Blue"
; gsn resources:
mpres@gsnStringFontHeightF = 0.015
mpres@gsnLeftString = str1 ; add left string
mpres@gsnRightString = str2 ; add right string
; Title resources:
mpres@tiMainString = title
mpres@tiMainOffsetYF = 0.0 ; Move the title down.
mpres@tiMainFontHeightF = 0.015
return(mpres)
end
begin ; main script
system("rm -f "+plotname+"tmp")
system("sed -n '/^FM/,$p' "+obfile+" > "+plotname+"tmp")
obfile = plotname+"tmp"
data = asciiread(obfile, -1, "string") ; -1 means read all rows.
cdata = stringtochar(data)
nline = dimsizes(cdata(:,0))
nsynop = 0
nmetar = 0
nship = 0
nbuoy = 0
nbogus = 0
nsound = 0
nairep = 0
npilot = 0
nsatem = 0
nsatob = 0
ngpspw = 0
ngpsrf = 0
nqscat = 0
nprofl = 0
nairsr = 0
nradar = 0
; find out the number of each observation type
do n = 0, nline-1 ; loop1 of parsing the file line by line
if ( cdata(n,0:5) .eq. "FM-12 " ) then
nsynop = nsynop + 1
end if
if ( cdata(n,0:5) .eq. "FM-15 " .or. cdata(n,0:5) .eq. "FM-16 ") then
nmetar = nmetar + 1
end if
if ( cdata(n,0:5) .eq. "FM-13 " ) then
nship = nship + 1
end if
if ( cdata(n,0:5) .eq. "FM-18 " .or. cdata(n,0:5) .eq. "FM-19 " ) then
nbuoy = nbuoy + 1
end if
if ( cdata(n,0:5) .eq. "FM-135" ) then
nbogus = nbogus + 1
end if
if ( cdata(n,0:5) .eq. "FM-35 " .or. cdata(n,0:5) .eq. "FM-36 ") then
nsound = nsound + 1
end if
if ( cdata(n,0:5) .eq. "FM-42 " .or. cdata(n,0:5) .eq. "FM-97 ") then
nairep = nairep + 1
end if
if ( cdata(n,0:5) .eq. "FM-32 " ) then
npilot = npilot + 1
end if
if ( cdata(n,0:5) .eq. "FM-86 " ) then
nsatem = nsatem + 1
end if
if ( cdata(n,0:5) .eq. "FM-88 " ) then
nsatob = nsatob + 1
end if
if ( cdata(n,0:5) .eq. "FM-111" ) then
ngpspw = ngpspw + 1
end if
if ( cdata(n,0:5) .eq. "FM-116" ) then
ngpsrf = ngpsrf + 1
end if
if ( cdata(n,0:5) .eq. "FM-281" ) then
nqscat = nqscat + 1
end if
if ( cdata(n,0:5) .eq. "FM-132" ) then
nprofl = nprofl + 1
end if
if ( cdata(n,0:5) .eq. "FM-133" ) then
nairsr = nairsr + 1
end if
if ( cdata(n,0:5) .eq. "FM-128" ) then
nradar = nradar + 1
end if
end do ; end loop1 of parsing the file line by line
if ( nsynop .gt. 0 )then
synop_lat = new(nsynop, float)
synop_lon = new(nsynop, float)
end if
if ( nmetar .gt. 0 )then
metar_lat = new(nmetar, float)
metar_lon = new(nmetar, float)
end if
if ( nship .gt. 0 )then
ship_lat = new(nship, float)
ship_lon = new(nship, float)
end if
if ( nbuoy .gt. 0 )then
buoy_lat = new(nbuoy, float)
buoy_lon = new(nbuoy, float)
end if
if ( nbogus .gt. 0 )then
bogus_lat = new(nbogus, float)
bogus_lon = new(nbogus, float)
end if
if ( nsound .gt. 0 )then
sound_lat = new(nsound, float)
sound_lon = new(nsound, float)
end if
if ( nairep .gt. 0 )then
airep_lat = new(nairep, float)
airep_lon = new(nairep, float)
end if
if ( npilot .gt. 0 )then
pilot_lat = new(npilot, float)
pilot_lon = new(npilot, float)
end if
if ( nsatem .gt. 0 )then
satem_lat = new(nsatem, float)
satem_lon = new(nsatem, float)
end if
if ( nsatob .gt. 0 )then
satob_lat = new(nsatob, float)
satob_lon = new(nsatob, float)
end if
if ( ngpspw .gt. 0 )then
gpspw_lat = new(ngpspw, float)
gpspw_lon = new(ngpspw, float)
end if
if ( ngpsrf .gt. 0 )then
gpsrf_lat = new(ngpsrf, float)
gpsrf_lon = new(ngpsrf, float)
end if
if ( nqscat .gt. 0 )then
qscat_lat = new(nqscat, float)
qscat_lon = new(nqscat, float)
end if
if ( nprofl .gt. 0 )then
profl_lat = new(nprofl, float)
profl_lon = new(nprofl, float)
end if
if ( nairsr .gt. 0 )then
airsr_lat = new(nairsr, float)
airsr_lon = new(nairsr, float)
end if
if ( nradar .gt. 0 )then
radar_lat = new(nradar, float)
radar_lon = new(nradar, float)
end if
isynop = 0
imetar = 0
iship = 0
ibuoy = 0
ibogus = 0
isound = 0
iairep = 0
ipilot = 0
isatem = 0
isatob = 0
igpspw = 0
igpsrf = 0
iqscat = 0
iprofl = 0
iairsr = 0
iradar = 0
do n = 0, nline-1 ; loop2 of parsing the file line by line
if ( cdata(n,0:5) .eq. "FM-12 " ) then
synop_lat(isynop) = stringtofloat(charactertostring(cdata(n,80:91)))
synop_lon(isynop) = stringtofloat(charactertostring(cdata(n,103:114)))
isynop = isynop + 1
end if
if ( cdata(n,0:5) .eq. "FM-15 " .or. cdata(n,0:5) .eq. "FM-16 ") then
metar_lat(imetar) = stringtofloat(charactertostring(cdata(n,80:91)))
metar_lon(imetar) = stringtofloat(charactertostring(cdata(n,103:114)))
imetar = imetar + 1
end if
if ( cdata(n,0:5) .eq. "FM-13 " ) then
ship_lat(iship) = stringtofloat(charactertostring(cdata(n,80:91)))
ship_lon(iship) = stringtofloat(charactertostring(cdata(n,103:114)))
iship = iship + 1
end if
if ( cdata(n,0:5) .eq. "FM-18 " .or. cdata(n,0:5) .eq. "FM-19 " ) then
buoy_lat(ibuoy) = stringtofloat(charactertostring(cdata(n,80:91)))
buoy_lon(ibuoy) = stringtofloat(charactertostring(cdata(n,103:114)))
ibuoy = ibuoy + 1
end if
if ( cdata(n,0:5) .eq. "FM-135" ) then
bogus_lat(ibogus) = stringtofloat(charactertostring(cdata(n,80:91)))
bogus_lon(ibogus) = stringtofloat(charactertostring(cdata(n,103:114)))
ibogus = ibogus + 1
end if
if ( cdata(n,0:5) .eq. "FM-35 " .or. cdata(n,0:5) .eq. "FM-36 ") then
sound_lat(isound) = stringtofloat(charactertostring(cdata(n,80:91)))
sound_lon(isound) = stringtofloat(charactertostring(cdata(n,103:114)))
isound = isound + 1
end if
if ( cdata(n,0:5) .eq. "FM-42 " .or. cdata(n,0:5) .eq. "FM-97 ") then
airep_lat(iairep) = stringtofloat(charactertostring(cdata(n,80:91)))
airep_lon(iairep) = stringtofloat(charactertostring(cdata(n,103:114)))
iairep = iairep + 1
end if
if ( cdata(n,0:5) .eq. "FM-32 " ) then
pilot_lat(ipilot) = stringtofloat(charactertostring(cdata(n,80:91)))
pilot_lon(ipilot) = stringtofloat(charactertostring(cdata(n,103:114)))
ipilot = ipilot + 1
end if
if ( cdata(n,0:5) .eq. "FM-86 " ) then
satem_lat(isatem) = stringtofloat(charactertostring(cdata(n,80:91)))
satem_lon(isatem) = stringtofloat(charactertostring(cdata(n,103:114)))
isatem = isatem + 1
end if
if ( cdata(n,0:5) .eq. "FM-88 " ) then
satob_lat(isatob) = stringtofloat(charactertostring(cdata(n,80:91)))
satob_lon(isatob) = stringtofloat(charactertostring(cdata(n,103:114)))
isatob = isatob + 1
end if
if ( cdata(n,0:5) .eq. "FM-111" ) then
gpspw_lat(igpspw) = stringtofloat(charactertostring(cdata(n,80:91)))
gpspw_lon(igpspw) = stringtofloat(charactertostring(cdata(n,103:114)))
igpspw = igpspw + 1
end if
if ( cdata(n,0:5) .eq. "FM-116" ) then
gpsrf_lat(igpsrf) = stringtofloat(charactertostring(cdata(n,80:91)))
gpsrf_lon(igpsrf) = stringtofloat(charactertostring(cdata(n,103:114)))
igpsrf = igpsrf + 1
end if
if ( cdata(n,0:5) .eq. "FM-281" ) then
qscat_lat(iqscat) = stringtofloat(charactertostring(cdata(n,80:91)))
qscat_lon(iqscat) = stringtofloat(charactertostring(cdata(n,103:114)))
iqscat = iqscat + 1
end if
if ( cdata(n,0:5) .eq. "FM-132" ) then
profl_lat(iprofl) = stringtofloat(charactertostring(cdata(n,80:91)))
profl_lon(iprofl) = stringtofloat(charactertostring(cdata(n,103:114)))
iprofl = iprofl + 1
end if
if ( cdata(n,0:5) .eq. "FM-133" ) then
airsr_lat(iairsr) = stringtofloat(charactertostring(cdata(n,80:91)))
airsr_lon(iairsr) = stringtofloat(charactertostring(cdata(n,103:114)))
iairsr = iairsr + 1
end if
if ( cdata(n,0:5) .eq. "FM-128" ) then
radar_lat(iradar) = stringtofloat(charactertostring(cdata(n,36:47)))
radar_lon(iradar) = stringtofloat(charactertostring(cdata(n,50:61)))
iradar = iradar + 1
end if
end do ; end loop2 of parsing the file line by line
markersize1 = 0.006
wks = gsn_open_wks(out_type, plotname)
gsn_define_colormap(wks,colormap)
nc1 = NhlNewColor(wks,.8,.8,.8) ; add light gray to colormap
gsres = True
gsres@gsMarkerIndex = 16 ; Use filled dots for markers.
gsres@gsMarkerColor = "navy"
gsres@gsMarkerSizeF = markersize1 ; default 0.007
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if ( proc_synop .and. nsynop .gt. 0 ) then
mpres1 = setmpres("SYNOP",date,nsynop+"")
map1 = gsn_csm_map(wks,mpres1)
dummy1 = gsn_add_polymarker(wks,map1,synop_lon(:),synop_lat(:),gsres)
; draw(map1)
end if
if ( proc_metar .and. nmetar .gt. 0 ) then
mpres2 = setmpres("METAR",date,nmetar+"")
map2 = gsn_csm_map(wks,mpres2)
dummy2 = gsn_add_polymarker(wks,map2,metar_lon(:),metar_lat(:),gsres)
end if
if ( proc_ship .and. nship .gt. 0 ) then
mpres3 = setmpres("SHIP",date,nship+"")
map3 = gsn_csm_map(wks,mpres3)
dummy3 = gsn_add_polymarker(wks,map3,ship_lon(:),ship_lat(:),gsres)
end if
if ( proc_sound .and. nsound .gt. 0 ) then
mpres4 = setmpres("sound",date,nsound+"")
map4 = gsn_csm_map(wks,mpres4)
dummy4 = gsn_add_polymarker(wks,map4,sound_lon(:),sound_lat(:),gsres)
end if
if ( proc_satob .and. nsatob .gt. 0 ) then
mpres5 = setmpres("SATOB",date,nsatob+"")
map5 = gsn_csm_map(wks,mpres5)
dummy5 = gsn_add_polymarker(wks,map5,satob_lon(:),satob_lat(:),gsres)
end if
resP = True
resP@gsnMaximize = True
resP@txString = "Observation"
gsn_panel(wks,(/dummy1,dummy2,dummy3,dummy4,dummy5/),(/2,3/),resP)
frame(wks)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Dec 7 02:33:25 2011
This archive was generated by hypermail 2.1.8 : Fri Dec 09 2011 - 13:09:09 MST