[ncl-talk] wrf domain plots

Tue May 3 06:39:28 MDT 2016

Dear NCL help,I am trying to plot the first domain of wrf with topography (also demarcating ocean and land area) and other two domains boundry embedded over it with the below script but i am unable to draw the boundry and also to get ocean land demarcation. I look forward for the support from ncl help in this regard. I also tried http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/Examples/EXPERIMENTAL/wrf_overlay_doms.htm this link but it always gives segmentation fault error.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" 
 begin   m1 = addfile("/mnt/gpfs1/data/private/fgvy024/Build_WRF/WPS/met_em.d01.2015-11-29_06:00:00.nc", "r")  m2 = addfile("/mnt/gpfs1/data/private/fgvy024/Build_WRF/WPS/met_em.d02.2015-11-29_06:00:00.nc", "r")  m3 = addfile("/mnt/gpfs1/data/private/fgvy024/Build_WRF/WPS/met_em.d03.2015-11-29_06:00:00.nc", "r")  m1t = wrf_user_getvar(m1,"HGT_M",0)  m2t = wrf_user_getvar(m2,"HGT_M",0)     ; Finer resolution  m3t = wrf_user_getvar(m3,"HGT_M",0) ; Finest resolution; Read lat/lon coordinates for each dataset.  m1lat = wrf_user_getvar(m1,"XLAT",0)   m1lon = wrf_user_getvar(m1,"XLONG",0)     m2lat = wrf_user_getvar(m2,"XLAT",0)   m2lon = wrf_user_getvar(m2,"XLONG",0)  m3lat = wrf_user_getvar(m3,"XLAT",0)  m3lon = wrf_user_getvar(m3,"XLONG",0)  wks = gsn_open_wks("png", "ccl")  ; send graphics to PNG file
  res = True   res at InitTime = False                        ; Do not plot time or footers  res at Footer = False   res at gsnDraw         = False   ; Don't draw plot (will do later)  res at gsnFrame        = False   ; Don't advance framce  (will do later);y axis levels  res at tiYAxisString = "Latitude"  res at tiYAxisFontHeightF = 0.015  res at tiYAxisFontColor = "black"  res at tiYAxisSide = "left"  res at tiYAxisAngleF = 90; x asis label  res at tiXAxisString = "Longitude"  res at tiXAxisFontHeightF = 0.015  res at tiXAxisFontColor = "black"   pltres = True  ;pltres at PanelPlot  = True  mpres1 = True  mpres1 at mpGeophysicalLineColor = "Black"  mpres1 at mpNationalLineColor    = "Black"  mpres1 at mpGridLineColor        = "Black"  mpres1 at mpLimbLineColor        = "Black"  mpres1 at mpPerimLineColor       = "Black"  mpres1 at mpFillOn = False  mpres1 at mpGridAndLimbOn       = False                ; Turn off lat/lon lines  mpres1 at pmTickMarkDisplayMode = "Always"  mpres1 at mpProjection          = "LambertConformal"    mpres1 at mpDataBaseVersion     = "MediumRes"          ; Default is LowRes  mpres1 at mpOutlineDrawOrder    = "PostDraw"           ; Draw map outlines last         opts_r = res     opts_r at cnFillOn = True    contour = wrf_contour(m1,wks,m1t,opts_r)  plot = wrf_map_overlays(m1,wks,contour,pltres,mpres1) ; Plot field over map background  overlay(plot,contour)   lnres1 = True   lnres1 at gsLineThicknessF = 1.8   lnres1 at gsLineColor = "white"   lnres1 at gsLineLabelFontColor = "red"   lnres1 at sLineLabelFontHeightF = 0.0145    ; Add some boxes to the map, showing the two finer-resolution; map areas.  dims = dimsizes(m2t)  r2   = dimsizes(dims)  r21  = dims(r2-1)  r22  = dims(r2-2)
  xbox = (/m2lon(0,0),m2lon(0,r21-1),m2lon(r22-1,r21-1), \           m2lon(r22-1,0),m2lon(0,0)/)  ybox = (/m2lat(0,0),m2lat(0,r21-1),m2lat(r22-1,r21-1), \           m2lat(r22-1,0),m2lat(0,0)/)
  x_out = new(dimsizes(xbox),typeof(xbox))  y_out = new(dimsizes(ybox),typeof(ybox))
; Can't use gsn_polyline here, because will get curved box lines.  datatondc(plot, xbox, ybox, x_out, y_out)  gsn_polyline_ndc(wks, x_out, y_out, lnres1)
  dims = dimsizes(m3t)  r3   = dimsizes(dims)  r31  = dims(r3-1)  r32  = dims(r3-2)
  xbox = (/m3lon(0,0),m3lon(0,r31-1),m3lon(r32-1,r31-1), \           m3lon(r32-1,0),m3lon(0,0)/)  ybox = (/m3lat(0,0),m3lat(0,r31-1),m3lat(r32-1,r31-1), \           m3lat(r32-1,0),m3lat(0,0)/)
  datatondc(plot, xbox, ybox, x_out, y_out)  gsn_polyline_ndc(wks, x_out, y_out, lnres1) draw(plot) frame(wks) end  

Regards
Ramchandra 

From: rammetro at hotmail.com
To: ncl-talk at ucar.edu
Subject: help ncl script for extraction of station data
Date: Tue, 19 Apr 2016 16:59:55 +0200

Dear NCL group,I tried to extract hourly air temperature data from wrf out and used the following scripts. would you figure out if there is any errors in it because when i run it it stops with showing time information and then pressing q several times it can execute the script but it is too time consuming and difficult to work because i have to press q for the size of times in file. I look foward for kind response

    ; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------
    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"    load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
    ; --------------  BEGINING OF NCL SCRIPT ----------------beginDATADir = "/mnt/gpfs1/data/private/fgvy024/Build_WRF/WRFV3/run/"FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d03* ")numFILES = dimsizes(FILES)a = addfiles(FILES+".nc","r");times = wrf_user_list_times(a)     ; get times in the filentimes = 4320            ; number of times in the filetemp = new(ntimes,float) do it = 0, ntimes-1                 ;Loop for the time: time = it;************************************************ ;  - Select lon & lat of the point of interest - ;************************************************  ; - The function wrf_user_ll_to_ij find the nearest grid point for a specific lat and lon
Latitude = 27.959167  ;phakdingLongitude = 86.81278res = True      res at returnInt = True                       point = wrf_user_ll_to_ij(a,Longitude,Latitude,res) x = point(0)-1y = point(1)-1;tc2 = wrf_user_getvar(a,"HGT",it)     ; heigth extractiontc2 = wrf_user_getvar(a,"T2",it)       ; temp extractiontc2 = tc2 -273.16
temp(it) = tc2(y,x) end do
    npts=ntimes    fName = "pyramid_1km.txt"    data  = new( npts, "string")                print("  Time      temp       ")         do it = 0, ntimes-1 ;ntimes-1
        print (sprintf("%5.0f",it)    +" " \          +sprintf("%21.2f", temp(it)) +"  " )
    end do                     ; end of time loop
    do it = 0,ntimes-1 ;ntimes-1
       data (it)= (sprintf("%5.0f",it)    +" " \          +sprintf("%21.2f", temp(it)) +"  "  )
 end do                     ; end of time loop asciiwrite (fName , data)end
Regards
Ramchandra KarkiPhD Student,  Institute of Geography, University of HamburgBundesstrase 55, Hamburg Germany

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160503/4680a220/attachment.html 