Kamal,
In your original message you had a subscript error in this line:
x = f->prcp(0:30,{36.06:49.44},{-92.64:-67.35})
So when you get integer grid indices from the getind_latlon2d
function, just substitute these indices into the same line, in place
of the original lat and lon values. You are changing the method used
from "coordinate subscripting" to "direct subscripting" with integer
indices. Also remove the curly brackets, which are only used with
coordinate subscripting.
ij = getind_latlon2d (lat2d, lon2d, lat, lon)
i1 = ij(0,0) ; (i, j) indices for first corner grid point
j1 = ij(0,1)
i2 = ij(1,0) ; (i, j) indices for second corner grid point
j2 = ij(1,1)
x = f->prcp(0:30,i1:i2,j1:j2) ; time and space subset of data array
This is untested, and it is possible that I have reversed the sense or
position of i and j. Please print intermediate values and try to
debug this part yourself.
You may also need to reverse the order of the two requested corner
points, or use the corner points from the opposite diagonal, to avoid
flipping the data array along the X or Y dimensions. To avoid
flipping, ensure that i1 < i2 and j1 < j2 before subscripting. This
can be automated with the min and max functions, but try it first
using the direct approach shown here, to get a feel for how this
works.
To make plots, you will probably also need matching subsets of the 2-D
coordinates. These subsets are extracted with exactly the same
subscripts as for the main data array:
lat2d_subset = lat2d(i1:i2,j1:j2)
lon2d_subset = lon2d(i1:i2,j1:j2)
Please review the basic NCL subscripting methods on this page. If you
are doing a lot of work with arrays and 2-D coordinates, you really
need a solid feel for this. HTH:
http://www.ncl.ucar.edu/Document/Manuals/Getting_Started/basics.shtml#ArraySubscripting
--Dave
On Tue, Apr 2, 2013 at 5:11 PM, <mmkamal@uwaterloo.ca> wrote:
> Hi Dave,
>
> First part of your answer is clear to me but failed to understand the second
> part. You stated that "you would pick two corners that are diagonally
> opposite in the desired final result, and translate the corner coordinates
> to integer grid indices.". I know the opposite corner coordinates of my area
> of interest which are as follows: (36.06,-92.64) & (49.44, -67.35). Now I
> follow you given link:
>
> ===========================================================================
> f = addfile("/brosto2/s3/kamal/DAYMET/prcp.sontario.2002.nc","r")
>
> lat2d = f->lat ; read 2-D coordinates
> lon2d = f->lon
>
> printVarSummary(lat2d)
> printMinMax(lat2d, True)
> printMinMax(lon2d, True)
>
>
> lat = (/ 36.06 , 49.44 /) ; user specified coordinate pairs
> lon = (/ -92.64, -67.35 /)
> ; return 2d subscripts
> nm = getind_latlon2d (lat2d,lon2d, lat, lon)
>
> print(nm)
>
> do k=0,dimsizes(lat)-1
> n = nm(k,0)
> m = nm(k,1)
> print(lat2d(n,m)+" "+lon2d(n,m))
>
> end do
>
> =================================================================
> Which gave me the following result:
>
> Variable: lat2d
> Type: double
> Total Size: 36571392 bytes
> 4571424 values
> Number of Dimensions: 2
> Dimensions and sizes: [y | 1924] x [x | 2376]
> Coordinates:
> y: [1241500..-681500]
> x: [153500..2528500]
> Number Of Attributes: 4
> standard_name : latitude
> long_name : latitude coordinate
> units : degrees_north
> actual_range : ( 32.06709695951115, 54.12263315213487 )
> (0)
> (0) latitude coordinate: min=32.0671 max=54.1226
> (0)
> (0) longitude coordinate: min=-98.2268 max=-62.8484
>
>
> Variable: nm
> Type: integer
> Total Size: 16 bytes
> 4 values
> Number of Dimensions: 2
> Dimensions and sizes: [2] x [2]
> Coordinates:
> Number Of Attributes: 1
> long_name : indices closest to specified LAT/LON coordinate pairs
> (0,0) 1897
> (0,1) 483
> (1,0) 66
> (1,1) 2060
> (0) 36.0613 -92.6386
> (0) 49.4393 -67.3466
> ======================================================================
>
> Now, I am confused what to do next. I look forward to hearing from you.
>
>
> Thanks
> Kamal
>
>
>
> Quoting Dave Allured - NOAA Affiliate <dave.allured@noaa.gov>:
>
>> Kamal,
>>
>> Oh, now I see the real problem. Normal NCL coordinate subscripting is
>> for 1-D coordinates, not 2-D coordinates such as those in Daymet
>> files. Your existing code is trying to subscript on the 1-D
>> projection coordinates X and Y (in meters), rather than the 2-D
>> geographic coordinates lat and lon (in degrees).
>>
>> There is a lookup function for 2-D coordinates. This should do what
>> you want. You input lat, lon coordinate pairs, and get back integer
>> indices which are then used as direct integer subscripts without the
>> curly braces. Typically, you would pick two corners that are
>> diagonally opposite in the desired final result, and translate the
>> corner coordinates to integer grid indices. Take a look at the
>> example on this function page:
>>
>>
>> http://www.ncl.ucar.edu/Document/Functions/Contributed/getind_latlon2d.shtml
>>
>> --Dave
>>
>> On Tue, Apr 2, 2013 at 3:31 PM, <mmkamal@uwaterloo.ca> wrote:
>>>
>>> Hi Dave,
>>>
>>> Thanks for pointing out the mistake. Yes, I also noticed that mistake
>>> after
>>> posting the query in NCL talk. I have corrected the mistake later and
>>> defined as follows but still get the same error message:
>>>
>>> f = addfile("/brosto2/s3/kamal/DAYMET/prcp.sontario.2002.nc","r")
>>> x = f->prcp(0:30,{36.06:49.44},{-92.64:-67.35})
>>> Rainfall = dim_avg_n_Wrap(x,(/1,2/))
>>>
>>>
>>> ================================================================================
>>> This is the same DAYMET data set I have generated using your contributed
>>> script and the header is as follows:
>>>
>>> [kamal@bro95 DAYMET]$ /work/kamal/brown/netcdf-4.0.1/bin/ncdump -h
>>> prcp.sontario.2002.nc
>>> netcdf prcp.sontario.2002 {
>>> dimensions:
>>> time = UNLIMITED ; // (365 currently)
>>> ntiles = 96 ;
>>> y = 1924 ;
>>> x = 2376 ;
>>> nv = 2 ;
>>> variables:
>>> short lambert_conformal_conic ;
>>> lambert_conformal_conic:grid_mapping_name =
>>> "lambert_conformal_conic" ;
>>> lambert_conformal_conic:longitude_of_central_meridian =
>>> -100. ;
>>> lambert_conformal_conic:latitude_of_projection_origin =
>>> 42.5
>>> ;
>>> lambert_conformal_conic:false_easting = 0. ;
>>> lambert_conformal_conic:false_northing = 0. ;
>>> lambert_conformal_conic:standard_parallel = 25., 60. ;
>>> int tiles_included(ntiles) ;
>>> tiles_included:actual_range = 11383, 12477 ;
>>> double lat(y, x) ;
>>> lat:standard_name = "latitude" ;
>>> lat:long_name = "latitude coordinate" ;
>>> lat:units = "degrees_north" ;
>>> lat:actual_range = 32.0670969595111, 54.1226331521349 ;
>>> double y(y) ;
>>> y:standard_name = "projection_y_coordinate" ;
>>> y:long_name = "y coordinate of projection" ;
>>> y:units = "m" ;
>>> y:actual_range = -681500., 1241500. ;
>>> double x(x) ;
>>> x:standard_name = "projection_x_coordinate" ;
>>> x:long_name = "x coordinate of projection" ;
>>> x:units = "m" ;
>>> x:actual_range = 153500., 2528500. ;
>>> double lon(y, x) ;
>>> lon:standard_name = "longitude" ;
>>> lon:long_name = "longitude coordinate" ;
>>> lon:units = "degrees_east" ;
>>> lon:actual_range = -98.2267604163355, -62.8484092336098 ;
>>> short yearday(time) ;
>>> yearday:long_name = "yearday" ;
>>> yearday:valid_range = 1s, 365s ;
>>> double time(time) ;
>>> time:long_name = "time" ;
>>> time:calendar = "standard" ;
>>> time:units = "days since 1980-01-01 00:00:00 UTC" ;
>>> time:bounds = "time_bnds" ;
>>> time:actual_range = 8036.5, 8400.5 ;
>>> double time_bnds(time, nv) ;
>>> float prcp(time, y, x) ;
>>> prcp:actual_range = 0.f, 200.f ;
>>> prcp:cell_methods = "area: sum time: sum" ;
>>> prcp:grid_mapping = "lambert_conformal_conic" ;
>>> prcp:coordinates = "lat lon" ;
>>> prcp:valid_range = 0.f, 200.f ;
>>> prcp:missing_value = -9999.f ;
>>> prcp:units = "mm/day" ;
>>> prcp:long_name = "daily total precipitation" ;
>>> prcp:_FillValue = -9999.f ;
>>>
>>> // global attributes:
>>> :history = "Fri Feb 1 23:30:46 EST 2013: Tiles joined by
>>> join.daymet.1226.ncl" ;
>>> :start_year = 2002s ;
>>> :source = "Daymet Software Version 2.0" ;
>>> :Version_software = "Daymet Software Version 2.0" ;
>>> :Version_data = "Daymet Data Version 2.1" ;
>>> :Conventions = "CF-1.4" ;
>>> :citation = "Please see http://daymet.ornl.gov/ for
>>> current
>>> Daymet data citation information" ;
>>> :references = "Please see http://daymet.ornl.gov/ for
>>> current information on Daymet references" ;
>>> }
>>>
>>> ======================================================================
>>>
>>> When I used x = f->prcp(0:30,{36.06:49.44},:) ; across all longitude then
>>> it
>>> works fine but getting error message with x =
>>> f->prcp(0:30,{36.06:49.44},{-92.64:-67.35}). So, it seems to me that my
>>> declaration of longitude ranges is not appropriate. I look forward to
>>> hearing from you.
>>>
>>>
>>>
>>> Thanks
>>> Kamal
>>>
>>>
>>>
>>>
>>> Quoting Dave Allured - NOAA Affiliate <dave.allured@noaa.gov>:
>>>
>>>> Kamal,
>>>>
>>>> In your first message you showed the same line of code twice. But the
>>>> first copy was different from the second, in a very important way:
>>>>
>>>> x = f->prcp(0:30,{36.06:49.44},{-92.64:-67.35-67.35}) ; read first 31
>>>> time
>>>> steps
>>>>
>>>> In this line, please notice that the number -67.35 was pasted TWICE
>>>> inside the curly brackets. I do not think you intended to do this.
>>>> Notice that with the way the line reads, "-67.35-67.35" is an
>>>> expression, not a simple number, and it evaluates to -134.7. If this
>>>> is in your actual code, then this will easily explain "Subscript out
>>>> of range". HTH.
>>>>
>>>> --Dave
>>>>
>>>> On Tue, Apr 2, 2013 at 12:17 PM, <mmkamal@uwaterloo.ca> wrote:
>>>>>
>>>>>
>>>>> Hi,
>>>>>
>>>>> I am trying to apply the NCL's "dim_avg_n_Wrap" function on a data
>>>>> set that contain data points beyond my are of interest.
>>>>>
>>>>> ; bounds over data present (Input Data)
>>>>>
>>>>> MinLat = 32.06
>>>>> MaxLat = 54.12
>>>>> MinLon = -98.22
>>>>> MaxLon = -62.84
>>>>>
>>>>> ; my are of interest
>>>>>
>>>>> MinLat = 36.06
>>>>> MaxLat = 49.44
>>>>> MinLon = -92.64
>>>>> MaxLon = -67.35
>>>>>
>>>>> I can use the following code to take average across all lat/lon
>>>>>
>>>>> =============================================================
>>>>>
>>>>> 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
>>>>> ;************************************************
>>>>> ; Read the file
>>>>> ;************************************************
>>>>>
>>>>> f = addfile("/brosto2/s3/kamal/DAYMET/prcp.sontario.2002.nc","r")
>>>>>
>>>>> x = f->prcp(0:30,{36.06:49.44},{-92.64:-67.35-67.35}) ; read first 31
>>>>> time steps
>>>>>
>>>>> Rainfall = dim_avg_n_Wrap(x,(/1,2/))
>>>>>
>>>>> end
>>>>> ==========================================================
>>>>>
>>>>> But I am getting the following error message:
>>>>>
>>>>> fatal:Subscript out of range, error in subscript #2
>>>>> fatal:Execute: Error occurred at or near line 19 in file
>>>>> Prcp_Jan_multi_time_labels_3.ncl
>>>>> ===========================================================
>>>>>
>>>>> Where the line 19 is:
>>>>>
>>>>> x = f->prcp(0:30,{36.06:49.44},{-92.64:-67.35})
>>>>>
>>>>>
>>>>> I have figure out that the problem is the declaration of longitude
>>>>> value. Could anyone please help me how can I define the longitude
>>>>> values.
>>>>>
>>>>>
>>>>> Thanks
>>>>> Kamal
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Apr 3 12:57:42 2013
This archive was generated by hypermail 2.1.8 : Thu Apr 04 2013 - 21:06:41 MDT