[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Fwd: Re: 20040322: Any idea why gribtocdl does a "floating exception" on this file?



Rich,

The cdl file is attached, it only has the true lat/lon for the 1st and the
last point. As I said before the calculations had too many rounding errors
or something. If you really need all those coordinates, you could create
them using MATLAB and then enter them into the cdl. It's only a text file.
Using  gribtonc to decode grib.3, it created a correct netcdf file. Hope
this helps. Been in meetings for the last couple of days.

Robb...




On Fri, 2 Apr 2004, Rich Signell wrote:

> Robb,
>
> Yes, my calcs were done using doubles (Matlab uses
> doubles by default).
>
> -Rich
>
> --- Robb Kambic <address@hidden> wrote:
> > On Thu, 1 Apr 2004, Richard Signell wrote:
> >
> > > Robb,
> > >
> > > I can't believe I didn't notice this before, but
> > you are using
> > > pole_lat, pole_lon = (-32.5, 10)
> > >
> > > but for the rtll formula I gave you to produce the
> > right results
> > > I was using
> > >
> > > pole_lat, pole_lon = (57.5, 10)
> > > which yields true lon/lat domain LL and UR
> > corners:
> > >
> > > lon(1,1)=4.5337
> > > lat(1,1)=33.2988
> > >
> > > lon(272,234)=24.6762
> > > lat(272,234)=49.4020
> >
> > Rich,
> >
> > The lats/lons are almost correct now, but I have
> > some round-off error.
> >
> > The first one is correct
> > la1 =-24, lo1 = -5
> > lon_true[ 1, 1 ] =4.53367311069887, lat_true[ 1, 1 ]
> > =33.2987847224989
> >
> > if I plug in the following values, then I get the
> > same answers you did.
> > la2 =-7.063, lo2 = 9.563
> > lon_true[ 2, 2 ] =24.67621972076, lat_true[ 2, 2 ]
> > =49.401952455077
> >
> > But iterations has round-off
> > lon[ 234 ] =10.269655932, lat[ 272 ] =-6.240013916
> > lon_true[ 234, 272 ] =26.0203076954498, lat_true[
> > 234, 272 ]
> > =50.0461679747345
> >
> > Next, I going to try doubles instead of floats.
> >
> > Robb...
> >
> >
> > >
> > > See attached image.
> > >
> > > So I guess either you need to insert a:
> > >
> > > pole_lat=90-pole_lat;
> > >
> > > line before you use my rtll formula, or else DWD
> > have
> > > incorrectly specified the "pole_lat" as
> > 90-"pole_lat" in the grib file.
> > >
> > > -Rich
> > >
> > >
> > > Robb Kambic wrote:
> > > >
> > > > On Tue, 30 Mar 2004, Richard Signell wrote:
> > > >
> > > > > Robb Kambic wrote:
> > > > > >
> > > > > > Rich,
> > > > > >
> > > > > > I have the formulas needed but the results
> > don't seem correct to me. The
> > > > > > pole is lat/lon  (-32.5, 10).  So ( -24, -5
> > ) -> ( -32.5, 10 ) and
> > > > > > (-7,063, 9.563 ) -> ( -15.563, 24.563 )
> > Since the rotation angle is 0,
> > > > > > there is no rotation.    The lat/lon
> > increment is .065535004
> > > > > >
> > > > > > The problem is the calculated corner from
> > the starting point of
> > > > > > ( -32.5, 10 ) ends up at ( -17.23, 27.76 )
> > instead of the given one
> > > > > > ( -15.563, 24.563 ).  It's a couple of
> > degrees off.  I don't think it's a
> > > > > > rounding err, so I must be missing
> > something.
> > > > > >
> > > > >
> > > > >
> > > > > I think I see your problem.   In the rotated
> > pole coordinate, the grid
> > > > > is uniformly spaced in lat/lon, but in the
> > rotated pole coordinate,
> > > > > it *is not*.   You therefore need 2D arrays to
> > hold the true lon/lat
> > > > > values, and you need to compute each of these
> > values using the formula.
> > > >
> > > > Rich,
> > > >
> > > > OK, I duplicated the code and the results still
> > are not correct. ie
> > > >
> > > > should be lon= 24.563 and lat =  -15.563
> > > >
> > > > straight increment, had error in first
> > calculation
> > > > lon[ 234 ] = 25.269655932, lat[ 272 ] =
> > -14.7400139159998
> > > >
> > > > results from duplicated rtll
> > > > almd=44.4907691654722, aphd=-43.1948333734209
> > > >
> > > > How about you send the lat lon values thru
> > MatLab and see what values you
> > > > get for the ending coordinates.
> > > >
> > > > Robb..
> > > >
> > > > >
> > > > > In other words, you need to loop thru each
> > grid point(e.g.):
> > > > >
> > > > > for j=1:272
> > > > >   lat=lat(j);
> > > > >   for i=1:234
> > > > >     lon=lon(i)
> > > > >
> >
> [lon_true(i,j),lat_true(i,j)]=rtll(pole_lon,pole_lat,lon(i),lat(i));
> > > > >   end
> > > > > end
> > > > >
> > > > > and lon_true, lat_true will be 2D arrays:
> > > > >
> > > > > lon_true(lat,lon)
> > > > > lat_true(lat,lon)
> > > > >
> > > > >
> > > > > > On Mon, 29 Mar 2004, Richard Signell wrote:
> > > > > >
> > > > > > > Robb,
> > > > > > >
> > > > > > > If it's helpful, here is my MATLAB
> > function for deriving true lon,lat
> > > > > > > from the rotated pole lon,lat.
> > > > > > >
> > > > > > > function
> > [almd,aphd]=rtll(tlm0d,tph0d,tlmd,tphd);
> > > > > > >
> >
> %-------------------------------------------------------------------------
> > > > > > > % RTLL transforms rotated coordinates
> > (tlmd,tphd) into
> > > > > > > % ordinary geographic coordinates
> > (almd,aphd). i/o decimal degrees (DD)
> > > > > > > % tlm0d, tph0d: lon e lat of the center of
> > rotation ("North Pole") in
> > > > > > > DD.
> > > > > > > %
> > > > > > > dtr=pi/180.;
> > > > > > >
> > > > > > > ctph0=cos(tph0d*dtr);
> > > > > > > stph0=sin(tph0d*dtr);
> > > > > > >
> > > > > > > stph=sin(tphd*dtr);
> > > > > > > ctph=cos(tphd*dtr);
> > > > > > > ctlm=cos(tlmd*dtr);
> > > > > > > stlm=sin(tlmd*dtr);
> > > > > > >
> > > > > > > aph=asin(stph0.*ctph.*ctlm+ctph0.*stph);
> > > > > > > cph=cos(aph);
> > > > > > >
> > > > > > > almd=tlm0d+asin(stlm.*ctph./cph)/dtr;
> > > > > > > aphd=aph/dtr;
> > > > > > > % end
> > > > > > >
> > > > > > >
> > > > > > > -Rich
> > > > > > >
> > > > > > >
> > > > >
> > > > > --
> > > > > Dr. Richard P. Signell                 |
> > address@hidden
> > > > > NATO/SACLANT Undersea Research Centre  |
> > Tel: (+39) 0187 527 381
> > > > > Viale San Bartolomeo 400               |
> > Fax: (+39) 0187 527 331
> > > > > 19138 La Spezia, ITALY  --> From USA/CANADA,
> > use: APO AE 09613-5000
> > > > >
> > > > >
> > > > > ***PRIVILEGED AND CONFIDENTIAL***
> > > > > The information contained in this e-mail
> > message (including any attached files) is intended
> > for the use of the addressee(s) only and is
> > priviliged information. The information should
> > neither be posted to the Internet, nor published in
> > any other public domain, without the express
> > permission
> === message truncated ===
>
>
> =====
> Rich Signell
> work: address@hidden
> home: address@hidden
>

===============================================================================
Robb Kambic                                Unidata Program Center
Software Engineer III                      Univ. Corp for Atmospheric Research
address@hidden             WWW: http://www.unidata.ucar.edu/
===============================================================================
netcdf Replace_with_model_name{ 


dimensions:
        record = UNLIMITED ;   // (reference time, forecast time)
        lat = 272 ;            // latitude
        lon = 234 ;            // longitude
        bls = 1 ;              // depth below land surface
        time_len = 21 ;        // string length for datetime strings
        valtime_offset = 7 ;   // number of offset times
        nmodels = 1 ;          // number of models
        ngrids = 1 ;           // number of grids
        nav = 1 ;              // for navigation
        nav_len = 100 ;        // max string length for navigation strings

variables:

        double reftime(record) ;        // reference time of the model
               reftime:long_name = "reference time" ;
               reftime:units = "hours since 1992-1-1" ;

        double valtime(record) ;        // forecast time ("valid" time)
               valtime:long_name = "valid time" ;
               valtime:units = "hours since 1992-1-1" ;

        :record = "reftime, valtime" ;  // "dimension attribute" -- means
                                        // (reftime, valtime) uniquely
                                        // determine record

        char   datetime(record, time_len) ; // derived from reftime
               datetime:long_name = "reference date and time" ;
               // units YYYY-MM-DD hh:mm:ssZ  (ISO 8601)

        double valtime_offset(valtime_offset) ; // valtime - reftime
               valtime_offset:long_name = "hours from reference time" ;
               valtime_offset:units = "hours" ;

        char   forecasttime(record, time_len) ; // derived from valtime
               forecasttime:long_name = "forecast date and time" ;
               // units YYYY-MM-DD hh:mm:ssZ  (ISO 8601)

        float  bls(bls) ;
               bls:long_name = "depth below land surface" ;
               bls:units = "cm" ;


        // The following lat and lon coordinate variables are redundant,
        // since the navigation variables provide the necessary information.
        // The extra information is included here for human readability.

//        float  lat(lat) ;
//               lat:long_name = "latitude" ;
//               lat:units = "degrees_north" ;

//        float  lon(lon) ;
//               lon:long_name = "longitude" ;
//               lon:units = "degrees_east" ;

        long   model_id(nmodels) ;
               model_id:long_name = "generating process ID number" ;

        // navigation variables all use nav dimension

        char   nav_model(nav, nav_len) ;        // navigation parameterization
               nav_model:long_name = "navigation model name" ;

        int    grid_type_code(nav) ;
               grid_type_code:long_name = "GRIB-1 GDS data representation type" 
;

        char   grid_type(nav, nav_len) ;
               grid_type:long_name = "GRIB-1 grid type" ;

        char   grid_name(nav, nav_len) ;
               grid_name:long_name = "grid name" ;

        int    grid_center(nav) ;
               grid_center:long_name = "GRIB-1 originating center ID" ;

        int    grid_number(nav, ngrids) ;
               grid_number:long_name = "GRIB-1 catalogued grid numbers" ;
               grid_number:_FillValue = -9999 ;

        char   i_dim(nav, nav_len) ;
               i_dim:long_name = "longitude dimension name" ;

        char   j_dim(nav, nav_len) ;
               j_dim:long_name = "latitude dimension name" ;

        int    Ni(nav) ;
               Ni:long_name = "number of points along a latitude circle" ;

        int    Nj(nav) ;
               Nj:long_name = "number of points along a longitude circle" ;

        float  La1(nav) ;
               La1:long_name = "latitude of first grid point" ;
               La1:units = "degrees_north" ;

        float  Lo1(nav) ;
               Lo1:long_name = "longitude of first grid point" ;
               Lo1:units = "degrees_east" ;

        float  La2(nav) ;
               La2:long_name = "latitude of last grid point" ;
               La2:units = "degrees_north" ;

        float  Lo2(nav) ;
               Lo2:long_name = "longitude of last grid point" ;
               Lo2:units = "degrees_east" ;

        float  La1_true(nav) ;
               La1_true:long_name = "true latitude of first grid point" ;
               La1_true:units = "degrees_north" ;

        float  Lo1_true(nav) ;
               Lo1_true:long_name = "true longitude of first grid point" ;
               Lo1_true:units = "degrees_east" ;

        float  La2_true(nav) ;
               La2_true:long_name = "true latitude of last grid point" ;
               La2_true:units = "degrees_north" ;

        float  Lo2_true(nav) ;
               Lo2_true:long_name = "true longitude of last grid point" ;
               Lo2_true:units = "degrees_east" ;

        float  Di(nav) ;
               Di:long_name = "longitudinal direction increment" ;
               Di:units = "degrees" ;

        float  Dj(nav) ;
               Dj:long_name = "latitudinal direction increment" ;
               Dj:units = "degrees" ;

        float  RotAngle(nav) ;
               RotAngle:long_name = "Angle of rotation" ;
               RotAngle:units = "degrees" ;

        float  RotLat(nav) ;
               RotLat:long_name = "Lat of S. pole of rotation" ;
               RotLat:units = "degrees" ;

        float  RotLon(nav) ;
               RotLon:long_name = "Lon of S. pole of rotation" ;
               RotLon:units = "degrees" ;

        byte   ResCompFlag(nav) ;
               ResCompFlag:long_name = "resolution and component flags" ;

        // end of navigation variables

        float  T_soil_bls(record,bls,lat,lon) ;
               T_soil_bls:long_name = "Soil temperature at depth below land 
surface" ;
               T_soil_bls:GRIB_parameter_number = 85 ;
               T_soil_bls:GRIB_level_flag = 111 ;
               T_soil_bls:units = "degK" ;
               T_soil_bls:_FillValue = -9999.f ;
               T_soil_bls:navigation = "nav" ;


// global attributes
               :history = "2004-03-29 12:37:57 - created by gribtocdl" ; 
               :title = "Enter model definition here" ;
               :Conventions = "NUWG" ;
               :GRIB_reference = "Office Note 388 GRIB" ;
               :GRIB_URL = "http://www.nco.ncep.noaa.gov/pmb/docs/on388/"; ;
               :version = 0.0 ;

data:

 bls = 0.0 ;
 model_id = 131 ;
 valtime_offset = 3, 6, 9, 12, 15, 18, 21 ;


 // Navigation
 nav_model = "GRIB1" ;
 grid_type_code = 10 ;
 grid_type = "Rotated latitude/longitude" ;
 grid_name = " " ;
 grid_center = 200 ;
 grid_number = 255 ;
 i_dim = "lon" ;
 j_dim = "lat" ;
 Ni = 234 ;
 Nj = 272 ;
 La1 = -24.000000 ;
 Lo1 = -5.000000 ;
 La2 = -7.063000 ;
 Lo2 = 9.563000 ;
 Di = 65.535004 ;
 Dj = 65.535004 ;
 RotLat = -32.500000 ;
 RotLon = 10.000000 ;
 RotAngle = 0.000000 ;
 La1_true = 33.2988 ;
 Lo1_true = 4.5337 ;
 La2_true = 49.4020 ;
 Lo2_true = 24.6762 ;
 ResCompFlag = 0 ;

}