[ncl-talk] Overlay profile of topography

Appo derbetini appopson4 at gmail.com
Tue Apr 25 00:56:41 MDT 2017


Good morning Mary,
Sorry , i forgot to add ncl-users email to the reply.

I am attaching the modified version of the script were I'm trying to create
a panel plot for two transects at latitude sellat and sellat2.
Unfortunately, only one plot appears.
Really, I don't know what is wrong this script.

Your help will be appreciated.



begin

    in1  = addfile("uvwr.mean.JJAS.nc", "r")                         ; open
netcdf file
    input4  = addfile("atanas_topo.nc", "r")

   u00     =  in1->u
   v00     =  in1->v
   w00     =  in1->w
   rhum00  =  in1->r

   topo1 = input4->HT


   lat = in1->latitude
   lon = in1->longitude


   dv00    = uv2dvF (u00, v00)
   uvd00   = dv2uvF (dv00)
   ucomp_div00 = uvd00(0, :, : , :, :)
   copy_VarCoords(u00, ucomp_div00)


   sellat = 3.0
   sellat2 = -2.5
   latmin = -5.0
   latmax = 3.0
   lonmin = min(lon)
   lonmax = 20.0

   udiv00 = dim_avg_n_Wrap(ucomp_div00(:, :, {latmin:latmax},
{lonmin:lonmax}), 2)
   omega00 = dim_avg_n_Wrap(w00(:, :, {latmin:latmax}, {lonmin:lonmax}), 2)
   r00 = dim_avg_n_Wrap(rhum00(:, :, {latmin:latmax},{lonmin:lonmax}), 2)

   ;udiv00  = ucomp_div00(:, :, {sellat}, {lonmin:lonmax})
   ;omega00 = w00(:, :, {sellat}, {lonmin:lonmax})
   ;r00     = rhum00(:, :, {sellat},{lonmin:lonmax})


   elev = 1013.25*(1 - topo1*0.0065/288.15)^5.25145

   copy_VarCoords(topo1, elev)

;  topo = dim_avg_n_Wrap(elev({latmin:latmax},{lonmin:lonmax}), 0)

   topo  = elev({sellat}, {lonmin:lonmax})
   topo2 = elev({sellat2}, {lonmin:lonmax})

;---create plot
   wks   = gsn_open_wks ("eps", "Overlay_Rhum_udiv_topo")        ; open
workstation
    plot = new(2, graphic)
   res                     = True                  ; plot mods desired
   res at gsnDraw             = False    ; turn off draw
   res at gsnFrame            = False    ; turn off frame
   res at gsnMaximize         = True

   res at gsnRightString      = ""
   res at gsnLeftString       = ""

   cnres = res
   cnres at cnFillOn          = True               ; turn on color fill
   cnres at cnFillPalette     = "MPL_rainbow"
   cnres at lbOrientation     = "Vertical"
   cnres at pmLabelBarOrthogonalPosF = 0.08
   cnres at tiYAxisString      = "Pressure (hPa)"

   vcres = res
   vcres at vcRefMagnitudeF = 10.0                ; define vector ref mag
   vcres at vcRefLengthF    = 0.1              ; define length of vec ref
   vcres at vcGlyphStyle    = "CurlyVector"      ; turn on curley vectors
   vcres at vcMinDistanceF  = 0.04               ; thin out vectors
   vcres at vcMapDirection  = False
   vcres at vcLineArrowThicknessF   = 3.0
   vcres at vcVectorDrawOrder  = "Draw"        ; draw vectors last

   contour_plot_1 = gsn_csm_pres_hgt(wks, r00(0, :, {lonmin:lonmax}), cnres)
   vector_plot_1  = gsn_csm_vector(wks, udiv00(0, :, {lonmin:lonmax}),
-200.0*omega00(0, :, {lonmin:lonmax}), vcres )

   contour_plot_2 = gsn_csm_pres_hgt(wks, r00(0, :, {lonmin:lonmax}), cnres)
   vector_plot_2  = gsn_csm_vector(wks, udiv00(0, :, {lonmin:lonmax}),
-200.0*omega00(0, :, {lonmin:lonmax}), vcres )


;---Add topo field using a filled polygon.
   getvalues contour_plot_1
      "trYMinF" : ymin
      "trYMaxF" : ymax
   end getvalues


;---Add topo field using a filled polygon.
   getvalues contour_plot_2
      "trYMinF" : ymin2
      "trYMaxF" : ymax2
   end getvalues



;---Create new X,Y arrays that form a closed polygon.
   nlon  = dimsizes(topo&lon)
   xtopo = new(nlon+3,typeof(topo&lon))
   ytopo = new(nlon+3,typeof(topo))

   xtopo(0:nlon-1) = topo&lon
   ytopo(0:nlon-1) = topo
   xtopo(nlon)     = topo&lon(nlon-1)
   ytopo(nlon)     = ymax          ; Use actual Y max of contour plot
   xtopo(nlon+1)   = topo&lon(0)
   ytopo(nlon+1)   = ymax
   xtopo(nlon+2)   = topo&lon(0)   ; This last point closes
   ytopo(nlon+2)   = topo(0)       ; the polygon.


   nlon2  = dimsizes(topo2&lon)
   xtopo2 = new(nlon2+3,typeof(topo2&lon))
   ytopo2 = new(nlon2+3,typeof(topo2))

   xtopo2(0:nlon-1) = topo2&lon
   ytopo2(0:nlon-1) = topo2
   xtopo2(nlon)     = topo2&lon(nlon-1)
   ytopo2(nlon)     = ymax2         ; Use actual Y max of contour plot
   xtopo2(nlon+1)   = topo2&lon(0)
   ytopo2(nlon+1)   = ymax2
   xtopo2(nlon+2)   = topo2&lon(0)   ; This last point closes
   ytopo2(nlon+2)   = topo2(0)       ; the polygon.


;---Add the polygon to the contour plot.
   gnres = True
   gnres at gsFillColor = "gray20"

   id = gsn_add_polygon(wks, contour_plot_1, xtopo, ytopo, gnres)
   id2 = gsn_add_polygon(wks, contour_plot_2, xtopo2, ytopo2, gnres)
   overlay(contour_plot_1, vector_plot_1)
   overlay(contour_plot_2, vector_plot_2)


     draw(contour_plot_1)     ; This draws everything
     draw(contour_plot_2)     ; This draws everything

   plot(0) = contour_plot_1

   plot(1) = contour_plot_2

    gsn_panel(wks,plot, (/2, 1/), res)               ; now draw as one plot


end


2017-04-24 21:58 GMT+01:00 Mary Haley <haley at ucar.edu>:

> Hi,
>
> In the future, please email ncl-talk with follow-up questions so everybody
> can benefit from the answers.
>
> What is the exact error message you're getting?
>
> --Mary
>
>
> On Mon, Apr 24, 2017 at 4:08 AM, Appo derbetini <appopson4 at gmail.com>
> wrote:
>
>> Hi Mary,
>>
>> Since you send me the script that plot masked fields, I am working to do
>> it for many latitude in a panel plot.
>> Always i got many errors saying that i cannot overlay.
>> How to solve it?
>> Thank you in advance.
>>
>> Appo
>>
>> 2017-04-10 8:23 GMT+01:00 Appo derbetini <appopson4 at gmail.com>:
>>
>>> Thank you very much.
>>> It's look like what I'm trying to do.
>>> regards
>>>
>>> 2017-04-09 22:11 GMT+01:00 Mary Haley <haley at ucar.edu>:
>>>
>>>> I'm sorry I didn't get back to this sooner, but I was having trouble
>>>> trying to understand what you are doing in the script.
>>>>
>>>> First, I'm not sure it makes sense to take an average of topographical
>>>> data across a set of latitudes and plot that as a topo layer.
>>>>
>>>> But, that aside, I think you just want to plot the topo line as a
>>>> filled polygon instead of a separate XY plot.
>>>>
>>>> I'm not sure what I've attached is correct, but hopefully it gives you
>>>> an idea of what to do.
>>>>
>>>> Note that I'm using "sellat" to select a single latitude value, rather
>>>> than averaging across a range of latitudes.
>>>>
>>>> --Mary
>>>>
>>>>
>>>> On Sun, Apr 9, 2017 at 2:10 AM, Appo derbetini <appopson4 at gmail.com>
>>>> wrote:
>>>>
>>>>> Dear Mary,
>>>>>
>>>>> Despite applying what you suggest the problem remain the same
>>>>> Thank you very much
>>>>>
>>>>> 2017-04-06 8:53 GMT+01:00 Appo derbetini <appopson4 at gmail.com>:
>>>>>
>>>>>> Dear Mary,
>>>>>>
>>>>>> Datasets are uploaded as indicated.
>>>>>>
>>>>>> What I really trying to do is something like Mask Example 14.
>>>>>>
>>>>>> But I don't want to interpolate topography to grid of wind.
>>>>>>
>>>>>> Here attached figure produced by the script
>>>>>>
>>>>>> Thank you
>>>>>>
>>>>>>
>>>>>> 2017-04-05 15:36 GMT+01:00 Mary Haley <haley at ucar.edu>:
>>>>>>
>>>>>>> I don't think I've ever seen this error before.
>>>>>>>
>>>>>>> It's hard to debug this particular script without having the data so
>>>>>>> I can run it. But, I'm wondering if this line might be part of the problem:
>>>>>>>
>>>>>>>   plot  = gsn_csm_pres_hgt_vector(wks, r00(0, :, {lonmin:lonmax}),
>>>>>>> udiv00(0, :, {lonmin:lonmax}), -200.0*omega00(0, :, {lonmin:lonmax}), res )
>>>>>>>
>>>>>>> When you do an arithmetic operation on an array, like:
>>>>>>>
>>>>>>>   -200. * omega00
>>>>>>>
>>>>>>> you end up stripping all of the metadata off omega00.  What's
>>>>>>> getting passed into gsn_csm_pres_hgt has no metadata. This might be okay,
>>>>>>> however, because I think this routine might just use the metadata from
>>>>>>> "udiv00" instead.
>>>>>>>
>>>>>>> Still, it's worth trying this:
>>>>>>>
>>>>>>>  omega00_scale = omega00(0, :, {lonmin:lonmax})   ; subsets omega00
>>>>>>> *and* copies all metadata
>>>>>>>  omega00_scale = -200.0 * omega00_scale
>>>>>>>
>>>>>>>   plot  = gsn_csm_pres_hgt_vector(wks, r00(0, :, {lonmin:lonmax}),
>>>>>>> udiv00(0, :, {lonmin:lonmax}), omega00_scale, res )
>>>>>>>
>>>>>>> If you continue to have problems with this script, could you upload
>>>>>>> your data (if it's not too large) to our ftp:
>>>>>>>
>>>>>>> ftp ftp.cgd.ucar.edu
>>>>>>> anonymous
>>>>>>> <use your email address for the password>
>>>>>>> cd incoming
>>>>>>> put uvwr.mean.JJAS.nc
>>>>>>> put atanas_topo.nc
>>>>>>> quit
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>> --Mary
>>>>>>>
>>>>>>>
>>>>>>> On Tue, Apr 4, 2017 at 9:18 AM, Appo derbetini <appopson4 at gmail.com>
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Dear all,
>>>>>>>> I am trying to overlay topography profil on a profil of wind and
>>>>>>>> relative humidity.
>>>>>>>> I don't want to interpolate topography to show some details.
>>>>>>>>
>>>>>>>> Atached are script and datasets used for the plot.
>>>>>>>>
>>>>>>>> But I having errors.
>>>>>>>>
>>>>>>>> warning:LlDataPolygon: point 6.150009,0.000000 outside data domain
>>>>>>>> warning:LlDataPolygon: point 19.983337,0.000000 outside data domain
>>>>>>>> warning:LlDataLineTo: point 6.150009,0.000000 outside data domain
>>>>>>>> warning:LlDataLineTo: point 19.983337,0.000000 outside data domain
>>>>>>>>
>>>>>>>> Any help will be appreciated.
>>>>>>>>
>>>>>>>> Best regards.
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> 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/20170425/17b77e0f/attachment.html 


More information about the ncl-talk mailing list