[ncl-talk] Fwd: How to add shapefiles to nesting domain
Tao Lu
hakufu.asano at gmail.com
Wed Jan 27 01:21:20 MST 2016
I try to show nesting domain and some shape files at the same time.
So, I edit the plotgrids_new.ncl as showed below, but I could not get what
I wanted.
I found the reason is the command: draw(mp), But I do not know how to
handle this problem.
Please help me.
; Script display location of model domains
; Only works for ARW domains
; Only works for NCL versions 6.2 or later
; 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.2) then
print("You need NCL V6.2 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 = "x11"
; 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 at mpFillOn = True
mpres at mpFillColors =
(/"background","DeepSkyBlue","ForestGreen","DeepSkyBlue", "transparent"/)
mpres at mpDataBaseVersion = "Ncarg4_1"
mpres at mpGeophysicalLineColor = "Black"
mpres at mpGridLineColor = "Black"
mpres at mpLimbLineColor = "Black"
mpres at mpNationalLineColor = "Black"
mpres at mpPerimLineColor = "Black"
mpres at mpUSStateLineColor = "Black"
; mpres at mpOutlineBoundarySets = "AllBoundaries"
;mpres at mpGridSpacingF = 45
mpres at tiMainString = " WPS Domain Configuration "
lnres = True
lnres at gsLineThicknessF = 2.5
lnres at domLineColors = (/ "white", "Red" , "Red" , "Blue" /)
txres = True
txres at txFont = "helvetica-bold"
;txres at txJust = "BottomLeft"
txres at txJust = "TopLeft"
txres at txPerimOn = False
txres at txFontHeightF = 0.015
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Do not change anything between the ";;;;;" lines
maxdom = 21
nvar = 19
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 at _FillValue = -999.0
plotvar = wrf_wps_read_nml(filename)
mpres at max_dom = floattointeger(plotvar(0,0))
mpres at dx = plotvar(0,1)
mpres at dy = plotvar(0,2)
if (.not.ismissing(plotvar(0,3))) then
mpres at ref_lat = plotvar(0,3)
else
mpres at ref_lat = 0.0
end if
if (.not.ismissing(plotvar(0,4))) then
mpres at ref_lon = plotvar(0,4)
else
mpres at ref_lon = 0.0
end if
if (.not.ismissing(plotvar(0,5))) then
mpres at ref_x = plotvar(0,5)
end if
if (.not.ismissing(plotvar(0,6))) then
mpres at ref_y = plotvar(0,6)
end if
mpres at truelat1 = plotvar(0,7)
mpres at truelat2 = plotvar(0,8)
mpres at stand_lon = plotvar(0,9)
mproj_int = plotvar(0,10)
mpres at pole_lat = plotvar(0,11)
mpres at pole_lon = plotvar(0,12)
do i = 0,maxdom-1
parent_idn(i) = floattointeger(plotvar(i,13))
parent_grid_ration(i) = floattointeger(plotvar(i,14))
i_parent_startn(i) = floattointeger(plotvar(i,15))
j_parent_startn(i) = floattointeger(plotvar(i,16))
e_wen(i) = floattointeger(plotvar(i,17))
e_snn(i) = floattointeger(plotvar(i,18))
end do
if(mpres at max_dom .gt. 1) then
do i = 1,mpres at 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 at 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 at parent_id = parent_idn(0:mpres at max_dom-1)
mpres at parent_grid_ratio = parent_grid_ration(0:mpres at max_dom-1)
mpres at i_parent_start = i_parent_startn(0:mpres at max_dom-1)
mpres at j_parent_start = j_parent_startn(0:mpres at max_dom-1)
mpres at e_we = e_wen(0:mpres at max_dom-1)
mpres at e_sn = e_snn(0:mpres at max_dom-1)
if(mproj_int .eq. 1) then
mpres at map_proj = "lambert"
mpres at pole_lat = 0.0
mpres at pole_lon = 0.0
else if(mproj_int .eq. 2) then
mpres at map_proj = "mercator"
mpres at pole_lat = 0.0
mpres at pole_lon = 0.0
else if(mproj_int .eq. 3) then
mpres at map_proj = "polar"
mpres at pole_lat = 0.0
mpres at pole_lon = 0.0
else if(mproj_int .eq. 4) then
mpres at map_proj = "lat-lon"
end if
end if
end if
end if
; Deal with global wrf domains that don't have dx or dy
if (mpres at dx.lt.1e-10 .and. mpres at dx.lt.1e-10) then
mpres at dx = 360./(mpres at e_we(0) - 1)
mpres at dy = 180./(mpres at e_sn(0) - 1)
mpres at ref_lat = 0.0
mpres at ref_lon = 180.0
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 at gsMarkerColor = "White"
;pmres at gsMarkerIndex = 16
;pmres at gsMarkerSizeF = 0.01
;gsn_polymarker(wks,mp,-77.26,38.56,pmres)
; Test to add a shapefile
; I add the following contents
dir = "/home/taolu/WRF/CASES/kinugawa/shp/"
shp_kinu = dir + "kinugawa.shp"
shp_tone = dir + "tonegawa.shp"
preres = True
preres at gsLineColor = "gray25"
preres at gsLineThicknessF = 2.0
kinu_id = gsn_add_shapefile_polylines(wks,mp,shp_kinu,preres)
tone_id = gsn_add_shapefile_polylines(wks,mp,shp_tone,preres)
draw(mp);
frame(wks) ; lets frame the plot - do not delete
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160127/e83b56e2/attachment.html
More information about the ncl-talk
mailing list