[ncl-talk] How to add shapefiles to nesting domain
Alan Brammer
abrammer at albany.edu
Wed Jan 27 12:22:26 MST 2016
> So, I edit the plotgrids_new.ncl as showed below, but I could not get what I wanted.
What do you mean by this?
Is a plot produced?
You’re outputting to x11 so it will only be temporary. At the top of the script you might want to change type =“png” or “pdf” then you can send the output to the script and detail what you actually wanted.
mpres at gsnDraw = False ; may prevent the wrf_wps_dom() function from drawing allowing your to overlay your shape file polylines. I’m not familiar with wrf_wps_dom() though, and the documentation I found online was particularly lacking.
An example of why the plot isn’t what you want, will enable someone to probably either answer or forward this onto wrfhelp at ucar.edu
Alan
##############################
Alan Brammer,
Post-Doc Researcher
Department of Atmospheric and Environmental Sciences,
University at Albany, State University of New York, Albany, NY, 12222
abrammer at albany.edu
##############################
> On 27 Jan 2016, at 03:21, Tao Lu <hakufu.asano at gmail.com> wrote:
>
> 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
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160127/6027b658/attachment.html
More information about the ncl-talk
mailing list