[ncl-talk] Help with creating a script for precipitation stored in shapefiles

[Carlos J. Valle-Diaz]

Hello,

I've been running WRF simulations to observe 24hr precipitation in Puerto
Rico from (July 12, 2011 12UTC to July 13, 2011 12UTC). I want to compare
the model runs with the observed from the National Weather Service. I have
their hourly shapefiles, see link to access them:
http://www.srh.noaa.gov/ridge2/Precip/qpehourlyshape/2011/201107/20110712/).

I managed to generate the hourly precipitation tendencies from the model
run by using a script from the NCL website and editing it accordingly to my
case (see attachment). Now, I want to do the same with the data reported by
the National Weather Service. I used a script provided by the NCL website
to print general information of the shapefiles and observed that the
geometry in the data set are points rather than polylines and polygons.
>From what I recall, there isn't an example script in the website to plot
hourly precipitation from shapefiles. If someone could provide me a basic
script, advice, and/or insight to generate the plots, it would be great. I
can take it from there and edit it accordingly to my case.

At the moment, I managed to generate in the script a map plot where Puerto
Rico is located. All that is left is to put the precipitation from the
shapefiles (see attachment).

I'm relatively new in NCL. I've seen the webinars and practiced with some
easy scripts. I appreciate any help with this! Whatever help you can
provide, I can take it from there. Also, I'm currently at NCAR as a
visiting scientist, so if the person that wants to help would like to meet,
he or she is very welcomed!

Cheers,
-- 
Carlos J. Valle Diaz
Ph.D. Chemistry Student
University of Puerto Rico Rio Piedras Campus
Department of Graduate Chemistry and Environmental Sciences
cj.vallediaz at gmail.com
Tel. 787-764-0000 x88192
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150902/f5d5d071/attachment.html 
-------------- next part --------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
  shp_name = nws_precip_2011071212.shp          ; File Name
  wks  = gsn_open_wks("X11","Rain PR-NWS")      ; Workstation

  res                     = True                ; Plot Options
  res at gsnDraw             = False
  res at gsnFrame            = False
  res at gsnMaximize         = True
  res at mpProjection        = "Mercator"
  res at mpLimitMode         = "LatLon"
  res at mpMinLatF           = 16.30
  res at mpMaxLatF           = 20.0
  res at mpMinLonF           = -64.0
  res at mpMaxLonF           = -68.0
  res at pmTickMarkDisplayMode = "Always"
  res at tiMainString          = "July 12-13, 2011 12UTC-24hr Precip"

  plot = gsn_csm_map(wks,res)                   ; Create Map

  f = addfile(shp_name, "r")                    ; Open Shapefile

end
---
-------------- next part --------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

; Get data
  a = addfile("./wrfout_d03_2011-07-12_00:00:00.nc","r")

; Workstation
  wks = gsn_open_wks("x11","plt_PrecipTend")

; Resources
  res = True
  res at MainTitle = "Hourly Precipitation: July 12-13, 2011 12UTC-24hour"
  pltres = True
  mpres = True
  mpres at mpGeophysicalLineColor = "Black"
  mpres at mpNationalLineColor    = "Black"
  mpres at mpUSStateLineColor     = "Black"
  mpres at mpGridLineColor        = "Black"
  mpres at mpLimbLineColor        = "Black"
  mpres at mpPerimLineColor       = "Black"

; Variable: Total precipitation
    rain_exp = wrf_user_getvar(a,"RAINNC",-1)
    rain_con = wrf_user_getvar(a,"RAINC",-1)
    rain_tot = rain_exp + rain_con
    rain_tot = rain_tot /25.4
 rain_tot at description = "Total Precipitation"
      rain_tot at units = "in"

; What times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)               ; number of times in the file

  nst = 12
  nend= nst+24

  do it = nst,nend                       ; Time Loop

    print("Working on time: " + times(it) )
    res at TimeLabel = times(it)   ; Set Valid time to use on plots

    rain_tot_tend = rain_tot(it,:,:) - rain_tot(nst,:,:)
    rain_tot_tend at description = "Precipitation Tendency"

   ; Plotting options for Precipitation
     opts_r = res
     opts_r at UnitLabel            = "in"
     opts_r at cnLevelSelectionMode = "ExplicitLevels"
     opts_r at cnLevels             = (/ .01, .1, .25, .5, 0.75, 1.0, 1.5, \
                                     2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10/)
     opts_r at cnFillColors         =(/"White","Cyan1","Cyan3","Blue","Green","Chartreuse3", \
                                    "ForestGreen","Yellow1","Yellow3","Orange", \
                                    "Red","Red4","HotPink1","Purple","Gray","Black"/)
     opts_r at cnInfoLabelOn        = False
     opts_r at cnConstFLabelOn      = False
     opts_r at cnFillOn             = True
     opts_r at SubFieldTitle = "from " + times(nst) + " to " + times(it)

   ; Total Precipitation (color fill)
     contour_tot = wrf_contour(a,wks, rain_tot_tend(:,:), opts_r)

   ; Make Plots

     ; Make Plot of Total Precipitation
       plot = wrf_map_overlays(a,wks,contour_tot,pltres,mpres)

  end do        ; End of Time Loop

end