Re: getind_latlon2d and xy plot

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Tue Mar 11 2014 - 09:43:16 MDT

Stephan,

Do you mean a contour plot, rather than an XY plot?

I don't what your data looks like (is it 2D, 3D, 4D?) or anything about your rotated grid, so it's hard to make suggestions.

Dave showed you how to get the coordinate locations in your 2D lat/lon grid for the area of interest. Using these coordinate locations, you can subscript the lat2d/lon2d arrays, *and* your original data array.

For example, say your data array is called "T", and is 3-dimensional (say, time x lat x lon), where the lat,lon are represented by those 2D lat2d/lon2d coordinates. Let's say you want to create a contour plot of the first timestep (index 0).

You could then subscript and plot your arrays as follows (UNTESTED):

;---Lat/lon area of interest
  minlat = -25
  maxlat = -10
  minlon = -75
  maxlon = -55

;---Get i,j index pairs closest to the given lat,lon values
  ji = getind_latlon2d(lat2d,lon2d,(/minlat,maxlat/),(/minlon,maxlon/))

  ilat1 = ji(0,0) ; Not necessary to assign to new
  ilat2 = ji(1,0) ; variables, but "ilat1" is
  ilon1 = ji(0,1) ; easier to remember than
  ilon2 = ji(1,1) ; "ji(0,0)".

  nt = 0 ; first timestep
  t_coch = t(nt,ilat1:ilat2,ilon1:ilon2)
  t_coch@lat2d = lat2d(ilat1:ilat2,ilon1:ilon2) ; necessary for plotting
  t_coch@lon2d = lon2d(ilat1:ilat2,ilon1:ilon2)

  wks = gsn_open_wks("x11","contour") ; "ps", "pdf", "png"

  res = True

  res@gsnMaximize = True ; maximize plot in frame
  res@cnFillOn = True ; turn on contour fill
  res@cnLinesOn = False ; turn off contour lines
  res@cnLineLabelsOn = False ; turn off line labels

;---Use if you have a high-res grid that could be slow to plot.
; res@cnFillMode = "RasterFill"

  res@tiMainString = "This is a main title"

  res@gsnAddCyclic = False ; because you are plotting regional data

  res@mpMinLatF = minlat
  res@mpMaxLatF = maxlat
  res@mpMinLonF = minlon
  res@mpMaxLonF = maxlon
  res@mpCenterLonF = (res@mpMinLonF + res@mpMaxLonF) / 2.

  plot = gsn_csm_contour_map(wks,t_coch,res)

n Mar 9, 2014, at 8:24 AM, Stephan Herrmann <stephan.w.herrmann@t-online.de> wrote:

> Hello everybody,
>
> how can I plot values from getind_latlon2d in a xy plot. I want to plot the mean temperature 1961-1970 for Cochabamba (17°S and 66°W). See my email (December 15th, 2013) and email from David Brown (December 16th, 2013) below. I think the solution ist quite simple but I don't know it.
>
>
> Best regards
>
> Stephan
>
>
>
> -----Original-Nachricht-----
> Betreff: Re: [ncl-talk] Error message concerning latitude and longitude
> Datum: Mon, 16 Dec 2013 21:20:07 +0100
> Von: David Brown <dbrown@ucar.edu>
> An: "Stephan Herrmann" <stephan.w.herrmann@t-online.de>
>
>
> Hi Stephan,
>
> If you uncomment the lines:
> ;res@mpMinLatF = min(lat2d) ; map range to zoom in on
> ;res@mpMaxLatF = max(lat2d)
> ;res@mpMinLonF = min(lon2d)
> ;res@mpMaxLonF = max(lon2d)
>
> And then change them to read:
> res@mpMinLatF = -25
> res@mpMaxLatF = -10
> res@mpMinLonF = -75
> res@mpMaxLonF = -55
> and then change
>
> res@mpLimitMode = "corners"
> to
> res@mpLimitMode = "latlon"
>
> you should get the the region of the map you want. Once you have the 2d coordinates you can draw the map in any projection with any map limits, and if the plot area is visible given your setting of the resources, it will display correctly.
>
>
> For 2) you can use
> ji = getind_latlon2d(lat2d,lon2d,-17,-66)
> print(ji)
> print ("closest grid point in data: " + lat2d(ji(0,0),ji(0,1)) + " " + lon2d(ji(0,0),ji(0,1)))
>
> This will give you:
>
> Variable: ji
> Type: integer
> Total Size: 8 bytes
> 2 values
> Number of Dimensions: 2
> Dimensions and sizes: [1] x [2]
> Coordinates:
> Number Of Attributes: 1
> long_name : indices closest to specified LAT/LON coordinate pairs
> (0,0) 52
> (0,1) 29
> (0) closest grid point in data: -16.9156 -66.1956
>
> Hope this helps.
> -dave
>
>
> On Dec 15, 2013, at 3:37 PM, "Stephan Herrmann" <stephan.w.herrmann@t-online.de> wrote:
>
> > Hello David,
> >
> > I have removed the lines printMinMax(lon2d,0) and printMinMax(lat2d,0) and now it works. I have two further questions:
> >
> > 1. I want to reduce the displayed map section (new map section: 10°S - 25°S and 55°W - 75°W). How can I do this?
> > 2. I want to plot the temperature for Cochabamba (17°S and 66°W) in a xy plot. How can I get the longitude and latitude?
> >
> >
> > Best regards
> >
> > Stephan
> >
> >
> >
> > -----Original-Nachricht-----
> > Betreff: Re: [ncl-talk] Error message concerning latitude and longitude
> > Datum: Sat, 14 Dec 2013 19:36:28 +0100
> > Von: David Brown <dbrown@ucar.edu>
> > An: Stephan Herrmann <stephan.w.herrmann@t-online.de>
> >
> > Hi Stephan,
> > Does the script actually quit? The min/max line is intended output from the printMinMax routine that I put in the script while debugging it. You could certainly remove it.
> > I don't understand where the "lines 1-1/1 (END)" comes from, but I do see that it was in your original message as well. Is it possible that you are using a modified version of some of the loaded script files (gsn_code.ncl, etc.)?
> > Otherwise, it is a strange message that is not in the usual format of NCL error messages. I suppose it's possible that some developer's personal debug message has been checked into the code accidentally -- that has happened occasionally.
> > -dave
> >
> > On Dec 14, 2013, at 9:04 AM, "Stephan Herrmann" <stephan.w.herrmann@t-online.de> wrote:
> >
> >> Hello David,
> >>
> >> thank you very much for your help. The PNG file of the graphics looks great. I also think the plot looks correct. When I run your NCL script test.ncl I get the following error message:
> >>
> >> (0) min=-83.0979 max=-51.1712
> >> lines 1-1/1 (END)
> >>
> >> How can I fix this problem?
> >>
> >>
> >> Best regards
> >>
> >> Stephan
> >>
> >>
> >>
> >> -----Original-Nachricht-----
> >> Betreff: Re: [ncl-talk] Error message concerning latitude and longitude
> >> Datum: Sat, 14 Dec 2013 02:00:12 +0100
> >> Von: David Brown <dbrown@ucar.edu>
> >> An: Stephan Herrmann <stephan.w.herrmann@t-online.de>
> >>
> >>
> >> Hi Stephan,
> >> I would like to add to the earlier suggestions from our team. The coordinates in the file are 1d rotated coordinates generated by CDO (I think). They are relative to the rotated coordinate system. For plotting purposes NCL needs to know the coordinates in the standard unrotated system, but currently there is no public interface for converting the coordinates. Also there are different conventions for specifying the rotated coordinates. I am attaching a little C source code that performs the operation based on the conventions employed by CDO. It generates an ascii output file containing the lat and lon of each grid point in the unrotated system. I am also attaching the ascii output, an NCL script that uses it to create the 2D coordinate arrays that locate your data, as well as a PNG file of the graphics. I think the plot looks correct, but let us know if it is not. Hopefully we can build this functionality into NCL in the near future.
> >> -dave
> >> _______________________________________________
> >> ncl-talk mailing list
> >> List instructions, subscriber options, unsubscribe:
> >> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> >
> >
> >
> >
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Mar 11 09:43:29 2014

This archive was generated by hypermail 2.1.8 : Fri Mar 14 2014 - 15:08:52 MDT