[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?
- Subject: Re: Fwd: Re: 20040322: Any idea why gribtocdl does a "floating exception" on this file?
- Date: Thu, 8 Apr 2004 14:34:18 -0600 (MDT)
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 ;
}