Re: How to calculate dim_avg_n_Wrap over selected latitude and longitude

From: <mmkamal_at_nyahnyahspammersnyahnyah>
Date: Tue Apr 02 2013 - 17:11:24 MDT

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 Tue Apr 2 17:11:43 2013

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