Re: Verifying a NARR grid value

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Thu, 11 Sep 2008 15:12:50 -0600 (MDT)

On Thu, 11 Sep 2008, Jamie Lahowetz wrote:

> Again thanks for all the help.
> I get this error:
> fatal:syntax error: function where expects 3 arguments, got 1
>
> I do admit I'm a little confused on what is being plotted. What I would like
> is the grid number (1,1 ; 30,30 ; etc if that is the format they are in)
> that the grid box represents this way I can compare the grid number that my
> program is grabbing with the actual grid number. I'm not sure which one does
> this.

I think I understand now. By grid numbers, I assume you mean the index
values of the lat/lon grid, and not the value of vpt.

If so, and you want to plot every grid number, then here's one way to
do it. It's not the most efficient, but it might help you see what's
going on. The '(i+1)+","+(j+1)' that you see below is taking the
number represented by 'i+1', appending a comma, and appending the
number represented by 'j+1', so you get a string that looks like
"1,1", "30,30", etc.

FYI: index subscripting in NCL starts at 0. It sounds like you want
"1,1" to be the label for the first grid point, so I'm accounting for
this in the code below.

    nlat = dimsizes(lat2d(:,0)) ; Get dimension sizes
    nlon = dimsizes(lon2d(0,:)) ; of the lat/lon dimensions.

    do i=0,nlat-1
      do j=0,nlon-1
        text_string = (i+1)+","+(j+1) ; "1,1", "1,2", etc.
        gsn_text(wks,plot,text_string,lon2d(i,j),lat2d(i,j),textres)
      end do
    end do

If you only want to print the grid number at certain lon,lat values,
like at lon = -101.01, lat = 39.60, then you need an extra "if"
statement:

    lat = 39.60
    lon = -101.01
    eps = 1e-3 ; or some small value
    do i=0,nlat-1
      do j=0,nlon-1
        if( (lat2d(i,j)-eps).le.lat.and.lat.le.(lat2d(i,j)+eps).and.\
            (lon2d(i,j)-eps).le.lon.and.lon.le.(lon2d(i,j)+eps)) then
          text_string = (i+1)+","+(j+1)
          gsn_text(wks,plot,text_string,lon2d(i,j),lat2d(i,j),textres)
        end if
      end do
    end do

Note that "do" loops can be expensive if you are looping across
big numbers. There are better ways of doing the code above, but
the above may be easier to follow.

If you still have problems, then if you can send me your data file and
your latest script, I can probably fix it for you.

--Mary

>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>
> ;############################
>
> begin
>
> fili =
> "/home/jlahowet/THOR/test2/narr/u-v_comp/narr-a_221_20080522_2100_000.grb"
> f = addfile (fili, "r")
>
> ;############################
> ; Print
> ;############################
>
> names = getfilevarnames(f)
> print(names)
> atts = getfilevaratts(f,names(0))
> dims = getfilevardims(f,names(0))
> print(atts)
> print(dims)
>
> ;############################
> ; Variables
> ;############################
>
> upt = f->USTM_221_HTGY
> vpt = f->VSTM_221_HTGY
> lat2d = f->gridlat_221
> lon2d = f->gridlon_221
>
> upt_at_lat2d = lat2d
> upt_at_lon2d = lon2d
> vpt_at_lat2d = lat2d
> vpt_at_lon2d = lon2d
>
> ;############################
>
> wks = gsn_open_wks("x11","lcnative")
> ;wks = gsn_open_wks("ps","lcnative")
> gsn_define_colormap(wks,"gui_default")
>
> ;############################
> ; Resources
> ;############################
>
> res = True
> res_at_tiMainString = "Native Lambert Conformal"
>
> res_at_mpLimitMode = "Corners"
> res_at_mpLeftCornerLatF = 37.0
> res_at_mpLeftCornerLonF = -103.0
> res_at_mpRightCornerLatF = 40.0
> res_at_mpRightCornerLonF = -98.0
>
> res_at_mpProjection = "LAMBERTCONFORMAL"
> res_at_mpLambertParallel1F = lon2d_at_Latin1
> res_at_mpLambertParallel2F = lon2d_at_Latin2
> res_at_mpLambertMeridianF = lat2d_at_Lov
>
> ; print(lon2d_at_Latin1)
> ; print(lon2d_at_Latin2)
> ; print(lon2d_at_Lov)
>
> res_at_pmTickMarkDisplayMode = "Always"
> res_at_mpFillOn = False
> res_at_mpOutlineDrawOrder = "PostDraw"
> res_at_mpOutlineBoundarySets = "GeophysicalAndUSStates"
> res_at_mpGridAndLimbOn = True
>
> res_at_tfDoNDCOverlay = False
>
> res_at_cnFillOn = True
> res_at_cnLinesOn = False
> res_at_gsnSpreadColors = True
> res_at_gsnAddCyclic = False
> res_at_gsnFrame = False
> res_at_gsnDraw = False
>
> ;############################
>
> plot = gsn_csm_contour_map(wks,upt,res)
> ;plot = gsn_csm_contour_map(wks,vpt,res)
>
> ;############################
> ; Contour
> ;############################
>
> cnres = True
> cnres_at_tfDoNDCOverlay = True ; don't transform
> cnres_at_cnInfoLabelOn = False ; turn off info label
> cnres_at_cnLineLabelsOn = True ; turn off line labels
> cnres_at_tiMainString = "Contour on Contour: Native Grid" ; title
> cnres_at_gsnDraw = False ; don't draw yet
> cnres_at_gsnFrame = False ; don't advance frame
> yet
>
> ;############################
>
> ;plotcn = gsn_contour(wks,vpt,cnres)
> ;overlay(plot,plotcn)
>
> ;############################
> ; Plot marker
> ;############################
>
> lon = -101.01
> lat = 39.60
> eps = 1e-3 ; or some small value
>
> polyres = True
> polyres_at_gsMarkerIndex = 16
> polyres_at_gsMarkerSizeF = 7.0
> polyres_at_gsMarkerColor = "black"
> draw(plot)
> gsn_polymarker(wks,plot,lon,lat,polyres)
>
> ;############################
> ; lat/lon Grid
> ;############################
>
> gridres = True
> gridres_at_gsMarkerIndex = 2
> gridres_at_gsMarkerSizeF = 10.0
> gridres_at_gsMarkerColor = "white"
> gsn_polymarker(wks,plot,ndtooned(lon2d),ndtooned(lat2d),gridres)
>
> ;############################
> ; Text Resources
> ;############################
>
> textres = True
> textres_at_txFontHeightF = 0.01 ; May need to change this
> textres_at_txFontColor = "black"
>
> ;############################
> ; value at every grid point
> ;############################
>
> ;lat1d = ndtooned(lat2d)
> ;lon1d = ndtooned(lon2d)
> ;vpt_text = ndtooned(vpt) + "" ; Convert to string
> ;gsn_text(wks,plot,vpt_text,lon1d,lat1d,textres)
>
> ;############################
> ; value at the specific point
> ;############################
>
> lat1d = ndtooned(lat2d) ; Convert to 1D
> lon1d = ndtooned(lon2d)
> ij = where( ((lat1d-eps).le.lat.and.lat.le.(lat1d+eps)).and. \
> ((lon1d-eps).le.lon.and.lon.le.(lon1d+eps)) )
> vpt_text = ndtooned(vpt) + "" ; Convert to string
> if(.not.any(ismissing(ij))) then
> gsn_text(wks,plot,vpt_text(ij),lon1d(ij),lat1d(ij),textres)
> end if
>
>
>
>
>
>
>
>
>
> ;vpt_text = ndtooned(vpt) + "" ; Convert to string
> ;vpt_text = sprintf("%7.2f",ndtooned(vpt))
> ;gsn_text(wks,plot,vpt_text,ndtooned(lon2d),ndtooned(lat2d),textres)
> ;gsn_text(wks,plot,vpt_text,ndtooned(x),ndtooned(y),textres)
>
> frame(wks)
>
> end
>
> On Thu, Sep 11, 2008 at 12:32 PM, Mary Haley <haley_at_ucar.edu> wrote:
>
>>
>> On Thu, 11 Sep 2008, Jamie Lahowetz wrote:
>>
>> How about if I wanted to print the grid number that is located at the
>>> polymarker (x,y)?
>>> With the script below I get the error:
>>> fatal:Number of dimensions on right hand side do not match number of
>>> dimension in left hand side
>>> fatal:Execute: Error occurred at or near line 4335 in file
>>> $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl
>>> fatal:Execute: Error occurred at or near line 135 in file
>>> 20080522-2100.ncl
>>>
>>>
>>>
>> You have:
>>
>> gsn_text(wks,plot,vpt_text,ndtooned(x),ndtooned(y),textres)
>>
>> which is effectively trying to plot a whole array of strings at a
>> single point. In order to use gsn_text, your x,y points must be the
>> same length as the array of text strings.
>>
>> If your desire is to give it some lat/lon point, and have it figure
>> out which value is at or close to that point, then you will need to
>> do something like this:
>>
>> lat1d = ndtooned(lat2d) ; Convert to 1D
>> lon1d = ndtooned(lon2d)
>> vpt_text = ndtooned(vpt) + "" ; Convert to string
>>
>> ; Lat/lon to locate
>> lon = -101.01
>> lat = 39.60
>> eps = 1e-3 ; or some small value
>>
>> ; Find where lon,lat is close or equal to lon1d,lat2d.
>> ij = where( ((lat1d-eps).le.lat.and.lat.le.(lat1d+eps)).and. \
>> ((lon1d-eps).le.lon.and.lon.le.(lon1d+eps)) )
>>
>> ; ij could be multiple indices.
>> if(.not.any(ismissing(ij))) then
>> gsn_text(wks,plot,vpt_text(ij),lon1d(ij),lat1d(ij),textres)
>> end if
>>
>> If you just want to plot the value of vpt at every lat2d/lon2d
>> location, then simply use:
>>
>> lat1d = ndtooned(lat2d)
>> lon1d = ndtooned(lon2d)
>> vpt_text = ndtooned(vpt) + "" ; Convert to string
>> gsn_text(wks,plot,vpt_text,lon1d,lat1d,textres)
>>
>> --Mary
>>
>>
>>
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>
>>> ;############################
>>>
>>> begin
>>>
>>> fili =
>>> "/home/jlahowet/THOR/test2/narr/u-v_comp/narr-a_221_20080522_2100_000.grb"
>>> f = addfile (fili, "r")
>>>
>>> ;############################
>>> ; Print
>>> ;############################
>>>
>>> names = getfilevarnames(f)
>>> print(names)
>>> atts = getfilevaratts(f,names(0))
>>> dims = getfilevardims(f,names(0))
>>> print(atts)
>>> print(dims)
>>>
>>> ;############################
>>> ; Variables
>>> ;############################
>>>
>>> upt = f->USTM_221_HTGY
>>> vpt = f->VSTM_221_HTGY
>>> lat2d = f->gridlat_221
>>> lon2d = f->gridlon_221
>>>
>>> upt_at_lat2d = lat2d
>>> upt_at_lon2d = lon2d
>>> vpt_at_lat2d = lat2d
>>> vpt_at_lon2d = lon2d
>>>
>>> ;############################
>>>
>>> ;wks = gsn_open_wks("x11","lcnative")
>>> wks = gsn_open_wks("ps","lcnative")
>>> gsn_define_colormap(wks,"gui_default")
>>>
>>> ;############################
>>> ; Resources
>>> ;############################
>>>
>>> res = True
>>> res_at_tiMainString = "Native Lambert Conformal"
>>>
>>> res_at_mpLimitMode = "Corners"
>>> res_at_mpLeftCornerLatF = 37.0
>>> res_at_mpLeftCornerLonF = -103.0
>>> res_at_mpRightCornerLatF = 40.0
>>> res_at_mpRightCornerLonF = -98.0
>>>
>>> res_at_mpProjection = "LAMBERTCONFORMAL"
>>> res_at_mpLambertParallel1F = lon2d_at_Latin1
>>> res_at_mpLambertParallel2F = lon2d_at_Latin2
>>> res_at_mpLambertMeridianF = lat2d_at_Lov
>>>
>>> ; print(lon2d_at_Latin1)
>>> ; print(lon2d_at_Latin2)
>>> ; print(lon2d_at_Lov)
>>>
>>> res_at_pmTickMarkDisplayMode = "Always"
>>> res_at_mpFillOn = False
>>> res_at_mpOutlineDrawOrder = "PostDraw"
>>> res_at_mpOutlineBoundarySets = "GeophysicalAndUSStates"
>>> res_at_mpGridAndLimbOn = True
>>>
>>> res_at_tfDoNDCOverlay = False
>>>
>>> res_at_cnFillOn = True
>>> res_at_cnLinesOn = False
>>> res_at_gsnSpreadColors = True
>>> res_at_gsnAddCyclic = False
>>> res_at_gsnFrame = False
>>> res_at_gsnDraw = False
>>>
>>> ;############################
>>>
>>> plot = gsn_csm_contour_map(wks,upt,res)
>>> ;plot = gsn_csm_contour_map(wks,vpt,res)
>>>
>>> ;############################
>>> ; Contour
>>> ;############################
>>>
>>> cnres = True
>>> cnres_at_tfDoNDCOverlay = True ; don't transform
>>> cnres_at_cnInfoLabelOn = False ; turn off info label
>>> cnres_at_cnLineLabelsOn = True ; turn off line labels
>>> cnres_at_tiMainString = "Contour on Contour: Native Grid" ;
>>> title
>>> cnres_at_gsnDraw = False ; don't draw yet
>>> cnres_at_gsnFrame = False ; don't advance frame
>>> yet
>>>
>>> ;############################
>>>
>>> ;plotcn = gsn_contour(wks,vpt,cnres)
>>> ;overlay(plot,plotcn)
>>>
>>> ;############################
>>> ; Plot marker
>>> ;############################
>>>
>>> x = -101.01
>>> y = 39.60
>>>
>>> polyres = True
>>> polyres_at_gsMarkerIndex = 16
>>> polyres_at_gsMarkerSizeF = 7.0
>>> polyres_at_gsMarkerColor = "black"
>>> draw(plot)
>>> gsn_polymarker(wks,plot,x,y,polyres)
>>>
>>> ;############################
>>> ; lat/lon Grid
>>> ;############################
>>>
>>> gridres = True
>>> gridres_at_gsMarkerIndex = 2
>>> gridres_at_gsMarkerSizeF = 10.0
>>> gridres_at_gsMarkerColor = "white"
>>> gsn_polymarker(wks,plot,ndtooned(lon2d),ndtooned(lat2d),gridres)
>>>
>>> ;############################
>>> ; Grid Value at point
>>> ;############################
>>>
>>> textres = True
>>> textres_at_txFontHeightF = 0.01 ; May need to change this
>>> textres_at_txFontColor = "black"
>>> vpt_text = ndtooned(vpt) + "" ; Convert to string
>>> ;vpt_text = sprintf("%7.2f",ndtooned(vpt))
>>> ;gsn_text(wks,plot,vpt_text,ndtooned(lon2d),ndtooned(lat2d),textres)
>>> gsn_text(wks,plot,vpt_text,ndtooned(x),ndtooned(y),textres)
>>>
>>> frame(wks)
>>>
>>> end
>>>
>>>
>>>
>>>
>>>
>>>
>>> On Thu, Sep 11, 2008 at 10:20 AM, Mary Haley <haley_at_ucar.edu> wrote:
>>>
>>>
>>>> Jamie,
>>>>
>>>> I believe the reason your marker is not showing up is because you are
>>>> using a longitude value that is outside the range of the map area you
>>>> are displaying.
>>>>
>>>> Your latitude/longitude area is defined by:
>>>>
>>>> res_at_mpLeftCornerLonF = -103.0
>>>> res_at_mpRightCornerLonF = -96.0
>>>> res_at_mpLeftCornerLatF = 36.0
>>>> res_at_mpRightCornerLatF = 43.0
>>>>
>>>> and you are trying to plot the point:
>>>>
>>>> x = 98
>>>> y = 56
>>>>
>>>> Longitude 98 is well outside this area, and latitude 56 is as well.
>>>>
>>>> By "grid" do you mean the lat/lon grid? If so, you could plot it using
>>>> the code you've set up below:
>>>>
>>>> polyres = True
>>>> polyres_at_gsMarkerIndex = 16
>>>> polyres_at_gsMarkerSizeF = 10.0
>>>> polyres_at_gsMarkerColor = "white" ; The (/ /) is not necessary
>>>> draw(plot)
>>>> gsn_polymarker(wks,plot,ndtooned(lon2d),ndtooned(lat2d),polyres)
>>>>
>>>> If you want to print the grid value at the lat/lon locations, then
>>>> you would use something like:
>>>>
>>>> textres = True
>>>> textres_at_txFontHeightF = 0.01 ; May need to change this
>>>> textres_at_txFontColor = "white"
>>>> draw(plot)
>>>> vpt_text = ndtooned(vpt) + "" ; Convert to string
>>>> gsn_text(wks,plot,vpt_text,ndtooned(lon2d),ndtooned(lat2d),textres)
>>>>
>>>> You can format the vpt strings using sprintf:
>>>>
>>>> vpt_text = sprintf("%5.2f",ndtooned(vpt))
>>>>
>>>> Of course, you may need to change the "5.2" format to accommodate
>>>> your data range.
>>>>
>>>> --Mary
>>>>
>>>>
>>>> On Wed, 10 Sep 2008, Jamie Lahowetz wrote:
>>>>
>>>> I'm looking for a way to verify that my algorithm is getting the right
>>>>
>>>>> grid
>>>>> points.
>>>>> Is there a way to display the grid and the number that that grid
>>>>> represents?
>>>>> Right now, The point I try to mark is not showing up and I used the same
>>>>> scheme as before, any ideas?
>>>>> Also, can I show the exact value that is represented by a grid position?
>>>>> Thanks for any help.
>>>>>
>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>>>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>>>>
>>>>> ;############################
>>>>>
>>>>> begin
>>>>>
>>>>> fili =
>>>>>
>>>>> "/home/jlahowet/THOR/test2/narr/u-v_comp/narr-a_221_20080522_2100_000.grb"
>>>>> f = addfile (fili, "r")
>>>>>
>>>>> ;############################
>>>>> ; Print
>>>>> ;############################
>>>>>
>>>>> names = getfilevarnames(f)
>>>>> print(names)
>>>>>
>>>>> atts = getfilevaratts(f,names(0))
>>>>> dims = getfilevardims(f,names(0))
>>>>> print(atts)
>>>>> print(dims)
>>>>>
>>>>> ;############################
>>>>> ; Variables
>>>>> ;############################
>>>>>
>>>>> upt = f->USTM_221_HTGY
>>>>> vpt = f->VSTM_221_HTGY
>>>>> lat2d = f->gridlat_221
>>>>> lon2d = f->gridlon_221
>>>>>
>>>>> upt_at_lat2d = lat2d
>>>>> upt_at_lon2d = lon2d
>>>>> vpt_at_lat2d = lat2d
>>>>> vpt_at_lon2d = lon2d
>>>>>
>>>>> ;############################
>>>>>
>>>>> wks = gsn_open_wks("x11","lcnative")
>>>>> gsn_define_colormap(wks,"gui_default")
>>>>>
>>>>> ;############################
>>>>> ; Resources
>>>>> ;############################
>>>>>
>>>>> res = True
>>>>> res_at_tiMainString = "Native Lambert Conformal"
>>>>>
>>>>> res_at_mpLimitMode = "Corners"
>>>>> res_at_mpLeftCornerLatF = 36.0
>>>>> res_at_mpLeftCornerLonF = -103.0
>>>>> res_at_mpRightCornerLatF = 43.0
>>>>> res_at_mpRightCornerLonF = -96.0
>>>>>
>>>>> res_at_mpProjection = "LAMBERTCONFORMAL"
>>>>> res_at_mpLambertParallel1F = lon2d_at_Latin1
>>>>> res_at_mpLambertParallel2F = lon2d_at_Latin2
>>>>> res_at_mpLambertMeridianF = lat2d_at_Lov
>>>>>
>>>>> ; print(lon2d_at_Latin1)
>>>>> ; print(lon2d_at_Latin2)
>>>>> ; print(lon2d_at_Lov)
>>>>>
>>>>> res_at_pmTickMarkDisplayMode = "Always"
>>>>> res_at_mpFillOn = False
>>>>> res_at_mpOutlineDrawOrder = "PostDraw"
>>>>> res_at_mpOutlineBoundarySets = "GeophysicalAndUSStates"
>>>>> res_at_mpGridAndLimbOn = True
>>>>>
>>>>> res_at_tfDoNDCOverlay = False
>>>>>
>>>>> res_at_cnFillOn = True
>>>>> res_at_cnLinesOn = False
>>>>> res_at_gsnSpreadColors = True
>>>>> res_at_gsnAddCyclic = False
>>>>> res_at_gsnFrame = False
>>>>> res_at_gsnDraw = False
>>>>>
>>>>> ;############################
>>>>>
>>>>> ;plot = gsn_csm_contour_map(wks,upt,res)
>>>>> plot = gsn_csm_contour_map(wks,vpt,res)
>>>>>
>>>>> ;############################
>>>>> ; Plot marker
>>>>> ;############################
>>>>>
>>>>> ;x = 39.60
>>>>> ;y = -101.01
>>>>> x = 98
>>>>> y = 56
>>>>> polyres = True
>>>>> polyres_at_gsMarkerIndex = 16
>>>>> polyres_at_gsMarkerSizeF = 10.0
>>>>> polyres_at_gsMarkerColor = (/"white"/)
>>>>> draw(plot)
>>>>> gsn_polymarker(wks,plot,x,y,polyres)
>>>>> frame(wks)
>>>>>
>>>>> end
>>>>>
>>>>>
>>>>>
>>>>> --
>>>>> Jamie Ryan Lahowetz
>>>>> University of Nebraska - Lincoln
>>>>> Graduate Student - Geosciences
>>>>> 402.304.0766
>>>>> jlahowe2_at_bigred.unl.edu
>>>>>
>>>>>
>>>>>
>>>
>>> --
>>> Jamie Ryan Lahowetz
>>> University of Nebraska - Lincoln
>>> Graduate Student - Geosciences
>>> 402.304.0766
>>> jlahowe2_at_bigred.unl.edu
>>>
>>>
>
>
> --
> Jamie Ryan Lahowetz
> University of Nebraska - Lincoln
> Graduate Student - Geosciences
> 402.304.0766
> jlahowe2_at_bigred.unl.edu
>
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Sep 11 2008 - 15:12:50 MDT

This archive was generated by hypermail 2.2.0 : Fri Sep 12 2008 - 21:42:04 MDT