Re: getind_latlon2d and xy plot

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Wed Mar 12 2014 - 09:19:07 MDT

Stephan,

Since you truly want an XY plot, then you don't need these two lines:

> T@lat2d = lat2d
> T@lon2d = lon2d

Then, in your call to gsn_csm_xy, you need to simply subscript T using the indices 52, 29:

> plot = gsn_csm_xy(wks,T&time,T(:,52,29),res)

If index 52 corresponds to longitude, and 29 to latitude, then you need to switch them:

> plot = gsn_csm_xy(wks,T&time,T(:,29,52),res)

Whenever you subscript a multi-dimensioned array with a single index value, then that dimension will no longer exist, and you end up with an array with one fewer dimensions.

So, the 3D array "T" becomes a 1D array when you do "T(:,29,52)" because you are subscripting *two* dimensions with a single index value.

The ":" means to keep all elements of that dimension, which is "time" in this case.

--Mary

On Mar 12, 2014, at 8:50 AM, Stephan Herrmann <stephan.w.herrmann@t-online.de> wrote:

> Hello Mary,
>
> thank you very much for your help! But I want to plot the annual temperature of Cochabamba (indices (52,29)) in a xy plot. My NCL code is as follows:
>
> 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"
>
> begin
>
> f = addfile("MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_yearmean.nc","r")
> T=f->var167
>
> latsize = getfilevardimsizes(f,"rlat")
> lonsize = getfilevardimsizes(f,"rlon")
> rotpole = f->rotated_pole ; need attributes for correct remap
>
> ll = asciiread("ll.out",(/latsize,lonsize,2/),"float") ; see email from David Brown
> (December 14th, 2013) below
>
> ;coordinates for Cochabamba
>
> lat2d = ll(52,29,0) ; insert indices from getind_latlon2d
> lon2d = ll(52,29,1) ; insert indices from getind_latlon2d
>
> T@lat2d = lat2d
> T@lon2d = lon2d
>
> wks=gsn_open_wks("png","xy_plot_temperature_Cochabamba")
> res = True
> res@tiYAxisString = "temperature (Kelvin)"
> res@tiXAxisString = "time (years)"
>
> res@trXMinF = 1961
> res@trXMaxF = 1990
> res@tiMainString = "Annual temperature of Cochabamba 1961-1990"
>
> plot = gsn_csm_xy(wks,T&time,annual temperature of Cochabamba (indices (52,29)),res)
>
> end
>
> What I have to write for annual temperature of Cochabamba (indices (52,29))?
>
> Here are the information about the variable "T":
>
> Variable: T
> Type: float
> Total Size: 556320 bytes
> 139080 values
> Number of Dimensions: 4
> Dimensions and sizes: [time | 30] x [height | 1] x [rlat | 76] x [rlon | 61]
> Coordinates:
> time: [19611231.95833333..19901231.95833333]
> height: [ 2.. 2]
> rlat: [-20.68..12.32]
> rlon: [157.56..183.96]
> Number Of Attributes: 1
> grid_mapping : rotated_pole
>
>
> Attachment:
>
> - file MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_yearmean.nc
> - file ll.out
>
> I Hope you can help me. I need this xy plot for my diploma thesis in Meteorology.
>
>
> Best regards
>
> Stephan
>
>
>
> -----Original-Nachricht-----
> Betreff: Re: [ncl-talk] getind_latlon2d and xy plot
> Datum: Tue, 11 Mar 2014 16:43:21 +0100
> Von: Mary Haley <haley@ucar.edu>
> An: Stephan Herrmann <stephan.w.herrmann@t-online.de>
>
>
> 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
> <MPIM-REMO_EH5OMA1Br3_HR_044_BOL_1961-1990_167_minus_1_hour_yearmean.nc><ll.out>_______________________________________________
> 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 Wed Mar 12 09:19:18 2014

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