[ncl-talk] How to add shapefiles to nesting domain
Alan Brammer
abrammer at albany.edu
Thu Jan 28 09:10:56 MST 2016
ok, I spent some time on this and possibly found a solution.
wrf_wps_dom() is calling draw() on the background map regardless of resources passed. It then overlays the domain plots using gsn_polyline_ndc(). This means that when you add your own polylines to mp and call draw() you are drawing over the domain boxes as they are in NDC coordinates and not attached to the actual plot.
If you read the lat/lon coordinates for your shape file (e.g. https://www.ncl.ucar.edu/Applications/Scripts/shapefiles_1.ncl <https://www.ncl.ucar.edu/Applications/Scripts/shapefiles_1.ncl>) . You could then use datatondc() https://www.ncl.ucar.edu/Document/Functions/Built-in/datatondc.shtml <https://www.ncl.ucar.edu/Document/Functions/Built-in/datatondc.shtml> and then add the outlines using gsn_polyline_ndc(wks, … … presres)
e.g.
n.b. lat and lon may not be stored as x/y in your shapefile you will need to investigate that yourself.
mp = wrf_wps_dom (wks,mpres,lnres,txres)
dir = "/home/taolu/WRF/CASES/kinugawa/shp/"
shp_kinu = dir + "kinugawa.shp"
f = addfile(dir + shp_kinu, "r”)
lon = f->x
lat = f->y
xndc = new(dimsizes(lat), float)
yndc = new(dimsizes(lat), float)
datatondc(mp, lon, lat, xndc, yndc)
gsn_polyline_ndc(wks, xndc, yndc, preres)
shp_tone = dir + "tonegawa.shp"
f = addfile(dir + shp_tone, "r”)
lon := f->x
lat := f->y
xndc := new(dimsizes(lat), float)
yndc := new(dimsizes(lat), float)
datatondc(mp, lon, lat, xndc, yndc)
gsn_polyline_ndc(wks, xndc, yndc, preres)
frame(wks)
;;; Don’t call draw(mp) at any time, this will cover up the domain boxes that were already added to the frame.
##############################
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 28 Jan 2016, at 10:29, Alan Brammer <abrammer at albany.edu> wrote:
>
> See reply below from Tao Lu. I don’t use the WRF functions enough to know whats going on in this cake.
>
> Unless obvious, always reply to ncl-talk not individual responders.
>
> Hopefully someone familiar with wrf_wps_dom() can comment on adding extra polylines to the output.
>
>
>
>>
>> Thank you very much for replying to me.
>> And I am sorry for my bad description.
>>
>> In the code I attached, only
>> ; 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);
>> is added by me.
>>
>> Run the original code, I could get the domains as showed below:
>> <image.png>
>>
>> But after I added my code, I could get this one:
>> <image.png>
>> The one I want is the combination of the two.
>>
>> I also tried type =“png” and mpres at gsnDraw = True, but also failed.
>> Do you have any clue about this?
>>
>> On Thu, Jan 28, 2016 at 4:22 AM, Alan Brammer <abrammer at albany.edu <mailto:abrammer at albany.edu>> wrote:
>>> 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 <mailto: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 <mailto:abrammer at albany.edu>
>> ##############################
>>
>>> On 27 Jan 2016, at 03:21, Tao Lu <hakufu.asano at gmail.com <mailto: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 <mailto:mpres at dx.lt>.1e-10 .and. mpres at dx.lt <mailto: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 <mailto:ncl-talk at ucar.edu>
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
>>
>>
>>
>>
>> --
>> ******************************************************
>> 盧 涛 (TAO LU)
>> 〒112-8551 東京都文京区春日1-13-27
>> 中央大学理工学研究科都市環境学専攻
>> 河川・水文研究室（山田正教授） 修士課程1年
>>
>> TEL: 03-3817-3406; Phone: 070-2188-7509
>> E-mail1: hakufu.asano at gmail.com <mailto:mail%3Amet.yamos at gmail.com>
>> E-mail2: lutao at civil.chuo-u.ac.jp <mailto:mail%3Ayamoto at civil.chuo-u.ac.jp>
>> *******************************************************
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160128/f6b852d6/attachment-0001.html
More information about the ncl-talk
mailing list