Re: How to calculate dim_avg_n_Wrap over selected latitude and longitude

From: Dave Allured - NOAA Affiliate <dave.allured_at_nyahnyahspammersnyahnyah>
Date: Tue Apr 02 2013 - 16:09:28 MDT

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 Tue Apr 2 16:09:39 2013

This archive was generated by hypermail 2.1.8 : Tue Apr 02 2013 - 21:23:48 MDT