- 积分
- 1619
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2019-3-13
- 最后登录
- 1970-1-1
|
30金钱
用WPS文件夹下的util里的plotgrids.ncl绘制domain区域图没有省市边界,ncl是6.4的,也尝试着用plotgrids_new.ncl画了,这个直接报错了,报错信息是fatal:Dimension sizes of left hand side and right hand side of assignment do not matchfatal:["Execute.c":8640]:Execute: Error occurred at or near line 75 in file plotgrids_new.ncl
plotgrids_new.ncl的程序如下 我基本没有改动过,我没啥要求,就想画个简单的嵌套区域图带着省市边界,这样我好判断我嵌套的区域是否合理,请各位帮忙看一下怎么处理合适,我是啥也不会
; Script display location of model domains
; Only works for ARW domains
; Only works for NCL versions later than 6.1
; Reads namelist file directly
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;
; Check the version of NCL
version = systemfunc("ncl -V")
if(version.lt.6.1) then
print("You need NCL V6.1 or later to run this script. Try running plotgrids.ncl. Stopping now...")
return
end if
; We generate plots, but what kind do we prefer?
type = "png"
; type = "pdf"
; type = "ps"
; type = "ncgm"
wks = gsn_open_wks(type,"wps_show_dom")
; read the following namelist file
filename = "namelist.wps"
; Set the colors to be used
colors = (/"white","black","White","ForestGreen","DeepSkyBlue","Red","Blue"/)
gsn_define_colormap(wks, colors)
; Set some map information ; line and text information
mpres = True
mpres@mpFillOn = True
mpres@mpFillColors = (/"background","DeepSkyBlue","ForestGreen","DeepSkyBlue", "transparent"/)
mpres@mpDataBaseVersion = "Ncarg4_1"
mpres@mpGeophysicalLineColor = "Black"
mpres@mpGridLineColor = "Black"
mpres@mpLimbLineColor = "Black"
mpres@mpNationalLineColor = "Black"
mpres@mpPerimLineColor = "Black"
mpres@mpUSStateLineColor = "Black"
; mpres@mpOutlineBoundarySets = "AllBoundaries"
;mpres@mpGridSpacingF = 45
mpres@tiMainString = " WPS Domain Configuration "
lnres = True
lnres@gsLineThicknessF = 2.5
lnres@domLineColors = (/ "white", "Red" , "Red" , "Blue" /)
txres = True
txres@txFont = "helvetica-bold"
;txres@txJust = "BottomLeft"
txres@txJust = "TopLeft"
txres@txPerimOn = False
txres@txFontHeightF = 0.015
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Do not change anything between the ";;;;;" lines
maxdom = 21
nvar = 17
parent_idn = new (maxdom,integer)
parent_grid_ration = new (maxdom,integer)
i_parent_startn = new (maxdom,integer)
j_parent_startn = new (maxdom,integer)
e_wen = new (maxdom,integer)
e_snn = new (maxdom,integer)
plotvar = new((/maxdom,nvar/),float)
plotvar@_FillValue = -999.0
plotvar = wrf_wps_read_nml(filename)
mpres@max_dom = floattointeger(plotvar(0,0))
mpres@dx = plotvar(0,1)
mpres@dy = plotvar(0,2)
mpres@ref_lat = plotvar(0,3)
mpres@ref_lon = plotvar(0,4)
mpres@truelat1 = plotvar(0,5)
mpres@truelat2 = plotvar(0,6)
mpres@stand_lon = plotvar(0,7)
mproj_int = plotvar(0,8)
mpres@pole_lat = plotvar(0,9)
mpres@pole_lon = plotvar(0,10)
do i = 0,maxdom-1
parent_idn(i) = floattointeger(plotvar(i,11))
parent_grid_ration(i) = floattointeger(plotvar(i,12))
i_parent_startn(i) = floattointeger(plotvar(i,13))
j_parent_startn(i) = floattointeger(plotvar(i,14))
e_wen(i) = floattointeger(plotvar(i,15))
e_snn(i) = floattointeger(plotvar(i,16))
end do
if(mpres@max_dom .gt. 1) then
do i = 1,mpres@max_dom-1
;Making sure edge is nested grid is at least 5 grid points from mother domain.
if(i_parent_startn(i) .lt. 5) then
print("Warning: Western edge of grid must be at least 5 grid points from mother domain!")
end if
if(j_parent_startn(i) .lt. 5) then
print("Warning: Southern edge of grid must be at least 5 grid points from mother domain!")
end if
pointwe = (e_wen(i)-1.)/parent_grid_ration(i)
pointsn = (e_snn(i)-1.)/parent_grid_ration(i)
gridwe = e_wen(parent_idn(i)-1)-(pointwe+i_parent_startn(i))
gridsn = e_snn(parent_idn(i)-1)-(pointsn+j_parent_startn(i))
if(gridwe .lt. 5) then
print("Warning: Eastern edge of grid must be at least 5 grid points from mother domain!")
end if
if(gridsn .lt. 5) then
print("Warning: Northern edge of grid must be at least 5 grid points from mother domain!")
end if
;Making sure nested grid is fully contained in mother domain.
gridsizewe = (((e_wen(parent_idn(i)-1)-4)-i_parent_startn(i))*parent_grid_ration(i))-(parent_grid_ration(i)-1)
gridsizesn = (((e_snn(parent_idn(i)-1)-4)-j_parent_startn(i))*parent_grid_ration(i))-(parent_grid_ration(i)-1)
if(gridwe .lt. 5) then
print("Warning: Inner nest (domain = " + (i+1) + ") is not fully contained in mother nest (domain = " + parent_idn(i) + ")!")
print("For the current setup of mother domain = " + parent_idn(i) + ", you can only have a nest of size " + gridsizewe + "X" + gridsizesn + ". Stopping Program!")
exit
end if
if(gridsn .lt. 5) then
print("Warning: Inner nest (domain = " + (i+1) + ") is not fully contained in mother nest (domain = " + parent_idn(i) + ")!")
print("For the current setup of mother domain = " + parent_idn(i) + ", you can only have a nest of size " + gridsizewe + "X" + gridsizesn + ". Stopping Program!")
exit
end if
;Making sure the nest ends of a mother grid domain point.
pointwetrunc = decimalPlaces(pointwe,0,False)
pointsntrunc = decimalPlaces(pointsn,0,False)
if((pointwe-pointwetrunc) .ne. 0.) then
nest_we_up = (ceil(pointwe)*parent_grid_ration(i))+1
nest_we_dn = (floor(pointwe)*parent_grid_ration(i))+1
print("Nest does not end on mother grid domain point. Try " + nest_we_dn + " or " + nest_we_up + ".")
end if
if((pointsn-pointsntrunc) .ne. 0.) then
nest_sn_up = (ceil(pointsn)*parent_grid_ration(i))+1
nest_sn_dn = (floor(pointsn)*parent_grid_ration(i))+1
print("Nest does not end on mother grid domain point. Try " + nest_sn_dn + " or " + nest_sn_up + ".")
end if
end do
end if
mpres@parent_id = parent_idn(0:mpres@max_dom-1)
mpres@parent_grid_ratio = parent_grid_ration(0:mpres@max_dom-1)
mpres@i_parent_start = i_parent_startn(0:mpres@max_dom-1)
mpres@j_parent_start = j_parent_startn(0:mpres@max_dom-1)
mpres@e_we = e_wen(0:mpres@max_dom-1)
mpres@e_sn = e_snn(0:mpres@max_dom-1)
if(mproj_int .eq. 1) then
mpres@map_proj = "lambert"
mpres@pole_lat = 0.0
mpres@pole_lon = 0.0
else if(mproj_int .eq. 2) then
mpres@map_proj = "mercator"
mpres@pole_lat = 0.0
mpres@pole_lon = 0.0
else if(mproj_int .eq. 3) then
mpres@map_proj = "polar"
mpres@pole_lat = 0.0
mpres@pole_lon = 0.0
else if(mproj_int .eq. 4) then
mpres@map_proj = "lat-lon"
end if
end if
end if
end if
mp = wrf_wps_dom (wks,mpres,lnres,txres)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Now you can add some information to the plot.
; Below is an example of adding a white dot over the DC location.
;pmres = True
;pmres@gsMarkerColor = "White"
;pmres@gsMarkerIndex = 16
;pmres@gsMarkerSizeF = 0.01
;gsn_polymarker(wks,mp,-77.26,38.56,pmres)
frame(wks) ; lets frame the plot - do not delete
end
|
最佳答案
查看完整内容
如果你懂NCL 可以直接在你的图上再覆盖一层边界线 边界线的来源:http://bbs.06climate.com/forum.php?mod=viewthread&tid=45252&extra=page%3D1
|