00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <sys/types.h>
00027 #include <sys/stat.h>
00028 #include <sys/mman.h>
00029 #include <unistd.h>
00030 #include <fcntl.h>
00031
00032 #include <Ebisuzaki.h>
00033 #include <GRIB.h>
00034
00035 #include <iostream>
00036 #include <sstream>
00037 #include <iomanip>
00038
00039
00040
00041
00042
00043 GRIB_TIME::GRIB_TIME( )
00044 {
00045 year = month = day = hour = minute = 0;
00046 p1 = p2 = p1p2 = naveraged = nmissing = 0;
00047 timeunit = GRIB_TIMEUNIT_UNKNOWN;
00048 timerange = GRIB_TIMERANGE_UNKNOWN;
00049 reftime = "2000-01-01 00:00:00 UTC";
00050 timestring = "Undefined time";
00051 }
00052
00053 void GRIB_TIME::strref( )
00054 {
00055 std::stringstream s;
00056 s << std::setw(4) << std::setfill('0') << year << "-"
00057 << std::setw(2) << std::setfill('0') << month << "-"
00058 << std::setw(2) << std::setfill('0') << day << " "
00059 << std::setw(2) << std::setfill('0') << hour << ":"
00060 << std::setw(2) << std::setfill('0') << minute << ":00 UTC";
00061 reftime = s.str( );
00062 return;
00063 }
00064
00065 void GRIB_TIME::strtime( )
00066 {
00067 std::stringstream s;
00068 char *unitstr = timeunitstr( );
00069 switch (timerange)
00070 {
00071 case GRIB_TIMERANGE_FORECAST_AT_REFTIME_PLUS_P1:
00072 if (p1 == 0)
00073 s << "Uninitialized analysis product or image for " << reftime;
00074 else
00075 s << "Forecast product valid for " << reftime
00076 << " + " << p1 << " " << unitstr;
00077 break;
00078 case GRIB_TIMERANGE_ANALYSIS_AT_REFTIME:
00079 s << "Initialized analysis product for " << reftime;
00080 break;
00081 case GRIB_TIMERANGE_VALID_IN_REFTIME_PLUS_P1_REFTIME_PLUS_P2:
00082 s << "Product valid between " << reftime
00083 << " + " << p1 << " " << unitstr << " and" << std::endl
00084 << " " << reftime
00085 << " + " << p2 << " " << unitstr;
00086 break;
00087 case GRIB_TIMERANGE_AVERAGE_IN_REFTIME_PLUS_P1_REFTIME_PLUS_P2:
00088 s << "Average from " << reftime
00089 << " + " << p1 << " " << unitstr << " and" << std::endl
00090 << " " << reftime << " + " << p2 << " " << unitstr;
00091 break;
00092 case GRIB_TIMERANGE_ACCUMULATED_INTERVAL_REFTIME_PLUS_P1_REFTIME_PLUS_P2:
00093 s << "Accumulated from " << reftime
00094 << " + " << p1 << " " << unitstr << " and" << std::endl
00095 << " " << reftime << " + " << p2 << " " << unitstr;
00096 break;
00097 case GRIB_TIMERANGE_DIFFERENCE_REFTIME_PLUS_P2_REFTIME_PLUS_P1:
00098 s << "Difference between product at " << reftime
00099 << " + " << p1 << " " << unitstr << " and" << std::endl
00100 << " product at "
00101 << reftime << " + " << p2 << " " << unitstr;
00102 break;
00103 case GRIB_TIMERANGE_AVERAGE_IN_REFTIME_MINUS_P1_REFTIME_MINUS_P2:
00104 s << "Average from " << reftime
00105 << " - " << p1 << " " << unitstr << " and" << std::endl
00106 << " " << reftime << " - " << p2 << " " << unitstr;
00107 break;
00108 case GRIB_TIMERANGE_AVERAGE_IN_REFTIME_MINUS_P1_REFTIME_PLUS_P2:
00109 s << "Average from " << reftime
00110 << " - " << p1 << " " << unitstr << " and" << std::endl
00111 << " " << reftime << " + " << p2 << " " << unitstr;
00112 break;
00113 case GRIB_TIMERANGE_VALID_AT_REFTIME_PLUS_P1P2:
00114 s << "Forecast product valid for " << reftime
00115 << " + " << p1p2 << " " << unitstr;
00116 break;
00117 case GRIB_TIMERANGE_CLIMATOLOGICAL_MEAN_OVER_MULTIPLE_YEARS_FOR_P2:
00118 s << "Climatological mean of " << naveraged << "/" << nmissing
00119 << " years from "
00120 << reftime << " for " << p2 << " " << unitstr;
00121 break;
00122 case GRIB_TIMERANGE_AVERAGE_OVER_FORECAST_OF_PERIOD_P1_REFTIME_PERIOD_P2:
00123 s << "Average of " << naveraged << "/" << nmissing;
00124 if (p1 == 0) s << " initialized analyses";
00125 else s << " forecasts with period " << p1 << " " << unitstr;
00126 s << std::endl << "starting at " << reftime << " at intervals of " << p2
00127 << " " << unitstr;
00128 break;
00129 case GRIB_TIMERANGE_ACCUMULATED_OVER_FORECAST_PERIOD_P1_REFTIME_PERIOD_P2:
00130 s << "Accumulated of " << naveraged << "/" << nmissing;
00131 if (p1 == 0) s << " initialized analyses";
00132 else s << " forecasts with period " << p1 << " " << unitstr;
00133 s << std::endl << "starting at " << reftime << " at intervals of " << p2
00134 << " " << unitstr;
00135 break;
00136 case GRIB_TIMERANGE_AVERAGE_OVER_FORECAST_OF_PERIOD_P1_AT_INTERVALS_P2:
00137 s << "Average of " << naveraged << "/" << nmissing
00138 << " forecasts at " << reftime
00139 << ", first with period " << p1 << " " << unitstr << std::endl
00140 << "and the other following at intervals of " << p2 << " " << unitstr;
00141 break;
00142 case GRIB_TIMERANGE_ACCUMULATION_OVER_FORECAST_PERIOD_P1_AT_INTERVALS_P2:
00143 s << "Accumulation of " << naveraged << "/" << nmissing
00144 << " forecasts at " << reftime
00145 << ", first with period " << p1 << " " << unitstr << std::endl
00146 << "and the other following at intervals of " << p2 << " " << unitstr;
00147 break;
00148 case GRIB_TIMERANGE_AVERAGE_OVER_FORECAST_FIRST_P1_OTHER_P2_REDUCED:
00149 s << "Average of " << naveraged << "/" << nmissing
00150 << " forecasts at " << reftime
00151 << ", first with" << std::endl << "period " << p1 << " " << unitstr
00152 << " and the other with period reduced by " << p2 << " " << unitstr;
00153 break;
00154 case GRIB_TIMERANGE_VARIANCE_OF_ANALYSES_WITH_REFERENCE_TIME_INTERVALS_P2:
00155 s << "Temporal variance or covariance of " << naveraged
00156 << "/" << nmissing
00157 << " initialized analyses," << std::endl << "starting at "
00158 << reftime << " with intervals of " << p2 << " " << unitstr;
00159 break;
00160 case GRIB_TIMERANGE_STDDEV_OF_FORECASTS_FIRST_P1_OTHER_P2_REDUCED:
00161 s << "Standard deviation of " << naveraged << "/" << nmissing
00162 << " forecasts at "
00163 << reftime << std::endl << "the first with period " << p1
00164 << " " << unitstr << " and the other following at "
00165 << p2 << " " << unitstr << " intervals";
00166 break;
00167 case GRIB_TIMERANGE_AVERAGE_OVER_ANALYSES_AT_INTERVALS_OF_P2:
00168 s << "Average of " << naveraged << "/" << nmissing
00169 << " uninitialized analyses "
00170 << "starting at " << reftime << "," << std::endl
00171 << "at intervals of " << p2 << " " << unitstr;
00172 break;
00173 case GRIB_TIMERANGE_ACCUMULATION_OVER_ANALYSES_AT_INTERVALS_OF_P2:
00174 s << "Accumulation of " << naveraged << "/" << nmissing
00175 << " uninitialized analyses "
00176 << "starting at " << reftime << "," << std::endl
00177 << "at intervals of " << p2 << " " << unitstr;
00178 break;
00179 case GRIB_TIMERANGE_STDDEV_OF_FORECASTS_RESPECT_TO_AVERAGE_OF_TENDENCY:
00180 s << "Standard deviation of " << naveraged << "/" << nmissing
00181 << " forecasts at "
00182 << reftime << ", with respect" << std::endl
00183 << "to time average of time tendency." << std::endl
00184 << "First forecast has period " << p1 << " " << unitstr
00185 << ", the other follow at intervals of " << p2 << " " << unitstr;
00186 break;
00187 case GRIB_TIMERANGE_AVERAGE_OF_DAILY_FORECAST_ACCUMULATIONS:
00188 s << "Average of " << naveraged << "/" << nmissing
00189 << " daily forecast accumulation"
00190 << " for " << reftime << std::endl << "from " << p1 << " " << unitstr
00191 << " to " << p2 << " " << unitstr << " with 1 day period.";
00192 break;
00193 case GRIB_TIMERANGE_AVERAGE_OF_SUCCESSIVE_FORECAST_ACCUMULATIONS:
00194 s << "Average of " << naveraged << "/" << nmissing
00195 << " successive forecast accumulation"
00196 << " for " << reftime << std::endl << "from " << p1 << " " << unitstr
00197 << " to " << p2 << " " << unitstr << " with " << (p2-p1) << " "
00198 << unitstr << " period";
00199 break;
00200 case GRIB_TIMERANGE_AVERAGE_OF_DAILY_FORECAST_AVERAGES:
00201 s << "Average of " << naveraged << "/" << nmissing
00202 << " daily forecast averages"
00203 << " for " << reftime << std::endl << "from " << p1 << " " << unitstr
00204 << " to " << p2 << " " << unitstr << " with 1 day period.";
00205 break;
00206 case GRIB_TIMERANGE_AVERAGE_OF_SUCCESSIVE_FORECAST_AVERAGES:
00207 s << "Average of " << naveraged << "/" << nmissing
00208 << " successive forecast averages"
00209 << " for " << reftime << std::endl << "from " << p1 << " " << unitstr
00210 << " to " << p2 << " " << unitstr << " with " << (p2-p1) << " "
00211 << unitstr << " period";
00212 default:
00213 std::cerr << "Unknown time range indicator: " << timerange << std::endl;
00214 break;
00215 }
00216 timestring = s.str( );
00217 return;
00218 }
00219
00220 void GRIB_TIME::Decode(unsigned char *pds)
00221 {
00222 year = PDS_Year4(pds);
00223 month = PDS_Month(pds);
00224 day = PDS_Day(pds);
00225 hour = PDS_Hour(pds);
00226 minute = PDS_Minute(pds);
00227
00228 timeunit = (t_enum_GRIB_TIMEUNIT) PDS_ForecastTimeUnit(pds);
00229 p1 = PDS_P1(pds);
00230 p2 = PDS_P2(pds);
00231 p1p2 = PDS_P1P2(pds);
00232 timerange = (t_enum_GRIB_TIMERANGE) PDS_TimeRange(pds);
00233 naveraged = PDS_NumAve(pds);
00234 nmissing = PDS_NumMissing(pds);
00235
00236 strref( );
00237 strtime( );
00238
00239 return;
00240 }
00241
00242 unsigned char *GRIB_TIME::Encode( )
00243 {
00244 static unsigned char time_encoded[13];
00245 time_encoded[0] = (year%100);
00246 time_encoded[1] = month;
00247 time_encoded[2] = day;
00248 time_encoded[3] = hour;
00249 time_encoded[4] = minute;
00250 time_encoded[5] = timeunit;
00251 time_encoded[6] = p1;
00252 time_encoded[7] = p2;
00253 time_encoded[8] = timerange;
00254 time_encoded[9] = (naveraged/256);
00255 time_encoded[10] = (naveraged%256);
00256 time_encoded[11] = nmissing;
00257 time_encoded[12] = (year/100)+1;
00258 return time_encoded;
00259 }
00260
00261 void GRIB_TIME::set(int year, int month, int day, int hour, int minute,
00262 t_enum_GRIB_TIMEUNIT tu, t_enum_GRIB_TIMERANGE tr,
00263 int time1, int time2, int nave, int nmiss)
00264 {
00265 set_referencetime(year, month, day, hour, minute);
00266 set_time(tu, tr, time1, time2, nave, nmiss);
00267 return;
00268 }
00269
00270 void GRIB_TIME::set_analysis( )
00271 {
00272 timeunit = GRIB_TIMEUNIT_HOUR;
00273 timerange = GRIB_TIMERANGE_ANALYSIS_AT_REFTIME;
00274 p1 = p2 = p1p2 = naveraged = nmissing = 0;
00275 strtime( );
00276 }
00277
00278 void GRIB_TIME::set_forecast_hour(int hforecast)
00279 {
00280 if (hforecast < 0)
00281 {
00282 std::cerr << "Negative or zero lenght forecast ????" << std::endl;
00283 throw;
00284 }
00285
00286 timeunit = GRIB_TIMEUNIT_HOUR;
00287 timerange = GRIB_TIMERANGE_FORECAST_AT_REFTIME_PLUS_P1;
00288 p1 = hforecast;
00289 p2 = p1p2 = 0;
00290 if (hforecast > 255)
00291 {
00292 timerange = GRIB_TIMERANGE_VALID_AT_REFTIME_PLUS_P1P2;
00293 p1p2 = hforecast;
00294 p1 = p1p2/256;
00295 p2 = p1p2%256;
00296 }
00297 naveraged = nmissing = 0;
00298 strtime( );
00299 }
00300
00301 void GRIB_TIME::set_referencetime(int year, int month, int day,
00302 int hour, int minute)
00303 {
00304 this->year = year;
00305 this->month = month;
00306 this->day = day;
00307 this->hour = hour;
00308 this->minute = minute;
00309 strref( );
00310 return;
00311 }
00312
00313 void GRIB_TIME::set_time(t_enum_GRIB_TIMEUNIT tu, t_enum_GRIB_TIMERANGE tr,
00314 int time)
00315 {
00316 this->set_time(tu, tr, time, 0, 0, 0);
00317 strtime( );
00318 return;
00319 }
00320
00321 void GRIB_TIME::set_time(t_enum_GRIB_TIMEUNIT tu, t_enum_GRIB_TIMERANGE tr,
00322 int time1, int time2)
00323 {
00324 this->set_time(tu, tr, time1, time2, 0, 0);
00325 strtime( );
00326 return;
00327 }
00328
00329 void GRIB_TIME::set_time(t_enum_GRIB_TIMEUNIT tu, t_enum_GRIB_TIMERANGE tr,
00330 int time1, int time2, int nave, int nmiss)
00331 {
00332 timeunit = tu;
00333 timerange = tr;
00334 switch (timerange)
00335 {
00336 case GRIB_TIMERANGE_FORECAST_AT_REFTIME_PLUS_P1:
00337 p1 = time1;
00338 break;
00339 case GRIB_TIMERANGE_ANALYSIS_AT_REFTIME:
00340 break;
00341 case GRIB_TIMERANGE_AVERAGE_IN_REFTIME_PLUS_P1_REFTIME_PLUS_P2:
00342 case GRIB_TIMERANGE_ACCUMULATED_INTERVAL_REFTIME_PLUS_P1_REFTIME_PLUS_P2:
00343 case GRIB_TIMERANGE_DIFFERENCE_REFTIME_PLUS_P2_REFTIME_PLUS_P1:
00344 case GRIB_TIMERANGE_AVERAGE_IN_REFTIME_MINUS_P1_REFTIME_MINUS_P2:
00345 case GRIB_TIMERANGE_AVERAGE_IN_REFTIME_MINUS_P1_REFTIME_PLUS_P2:
00346 p1 = time1;
00347 p2 = time2;
00348 break;
00349 case GRIB_TIMERANGE_VALID_AT_REFTIME_PLUS_P1P2:
00350 p1 = time1/256;
00351 p2 = time1%256;
00352 p1p2 = time1;
00353 break;
00354 case GRIB_TIMERANGE_CLIMATOLOGICAL_MEAN_OVER_MULTIPLE_YEARS_FOR_P2:
00355 case GRIB_TIMERANGE_AVERAGE_OVER_FORECAST_OF_PERIOD_P1_REFTIME_PERIOD_P2:
00356 case GRIB_TIMERANGE_ACCUMULATED_OVER_FORECAST_PERIOD_P1_REFTIME_PERIOD_P2:
00357 case GRIB_TIMERANGE_AVERAGE_OVER_FORECAST_OF_PERIOD_P1_AT_INTERVALS_P2:
00358 case GRIB_TIMERANGE_ACCUMULATION_OVER_FORECAST_PERIOD_P1_AT_INTERVALS_P2:
00359 case GRIB_TIMERANGE_AVERAGE_OVER_FORECAST_FIRST_P1_OTHER_P2_REDUCED:
00360 case GRIB_TIMERANGE_STDDEV_OF_FORECASTS_FIRST_P1_OTHER_P2_REDUCED:
00361 case GRIB_TIMERANGE_STDDEV_OF_FORECASTS_RESPECT_TO_AVERAGE_OF_TENDENCY:
00362 case GRIB_TIMERANGE_AVERAGE_OF_DAILY_FORECAST_ACCUMULATIONS:
00363 case GRIB_TIMERANGE_AVERAGE_OF_DAILY_FORECAST_AVERAGES:
00364 case GRIB_TIMERANGE_AVERAGE_OF_SUCCESSIVE_FORECAST_ACCUMULATIONS:
00365 case GRIB_TIMERANGE_AVERAGE_OF_SUCCESSIVE_FORECAST_AVERAGES:
00366 p1 = time1;
00367 p2 = time2;
00368 naveraged = nave;
00369 nmissing = nmiss;
00370 break;
00371 case GRIB_TIMERANGE_VARIANCE_OF_ANALYSES_WITH_REFERENCE_TIME_INTERVALS_P2:
00372 case GRIB_TIMERANGE_AVERAGE_OVER_ANALYSES_AT_INTERVALS_OF_P2:
00373 case GRIB_TIMERANGE_ACCUMULATION_OVER_ANALYSES_AT_INTERVALS_OF_P2:
00374 p1 = 0;
00375 p2 = time2;
00376 naveraged = nave;
00377 nmissing = nmiss;
00378 break;
00379 default:
00380 std::cerr << "Unknown time range indicator: " << timerange << std::endl;
00381 break;
00382 }
00383 strtime( );
00384 return;
00385 }
00386
00387 char *GRIB_TIME::timeunitstr( )
00388 {
00389 switch (timeunit)
00390 {
00391 case GRIB_TIMEUNIT_MINUTE:
00392 return "minute";
00393 break;
00394 case GRIB_TIMEUNIT_HOUR:
00395 return "hour";
00396 break;
00397 case GRIB_TIMEUNIT_DAY:
00398 return "day";
00399 break;
00400 case GRIB_TIMEUNIT_MONTH:
00401 return "month";
00402 break;
00403 case GRIB_TIMEUNIT_YEAR:
00404 return "year";
00405 break;
00406 case GRIB_TIMEUNIT_DECADE:
00407 return "decade";
00408 break;
00409 case GRIB_TIMEUNIT_NORMAL:
00410 return "normal";
00411 break;
00412 case GRIB_TIMEUNIT_CENTURY:
00413 return "century";
00414 break;
00415 case GRIB_TIMEUNIT_HOURS3:
00416 return "x3 hours";
00417 break;
00418 case GRIB_TIMEUNIT_HOURS6:
00419 return "x6 hours";
00420 break;
00421 case GRIB_TIMEUNIT_HOURS12:
00422 return "x12 hours";
00423 break;
00424 case GRIB_TIMEUNIT_SECOND:
00425 return "second";
00426 break;
00427 default:
00428 return "unknown unit";
00429 break;
00430 }
00431 }
00432
00433 time_t GRIB_TIME::ForecastSeconds( )
00434 {
00435 time_t rincr;
00436 int dtime = p1;
00437
00438 if (timerange == GRIB_TIMERANGE_VALID_AT_REFTIME_PLUS_P1P2)
00439 dtime = p1p2;
00440 if (timerange > GRIB_TIMERANGE_ANALYSIS_AT_REFTIME &&
00441 timerange < GRIB_TIMERANGE_AVERAGE_IN_REFTIME_MINUS_P1_REFTIME_MINUS_P2)
00442 dtime = p2;
00443
00444 if (timeunit == GRIB_TIMEUNIT_MINUTE)
00445 rincr = dtime * 60;
00446 else if (timeunit == GRIB_TIMEUNIT_HOUR)
00447 rincr = dtime * 3600;
00448 else if (timeunit == GRIB_TIMEUNIT_DAY)
00449 rincr = dtime * 86400;
00450 else if (timeunit == GRIB_TIMEUNIT_MONTH)
00451 rincr = dtime * 2592000;
00452 else if (timeunit == GRIB_TIMEUNIT_HOURS3)
00453 rincr = dtime * 10800;
00454 else if (timeunit == GRIB_TIMEUNIT_HOURS6)
00455 rincr = dtime * 21600;
00456 else if (timeunit == GRIB_TIMEUNIT_HOURS12)
00457 rincr = dtime * 43200;
00458 else if (timeunit == GRIB_TIMEUNIT_SECOND)
00459 rincr = dtime;
00460 else
00461 rincr = 0;
00462 return rincr;
00463 }
00464
00465 std::string GRIB_TIME::Reftime( ) { return reftime; }
00466 std::string GRIB_TIME::TimeString( ) { return timestring; }
00467
00468 std::string GRIB_TIME::Reftime2000( )
00469 {
00470 std::stringstream s;
00471 s << "2000-01-01 00:00:00 UTC";
00472 return s.str( );
00473 }
00474
00475 time_t GRIB_TIME::ForecastSeconds2000( )
00476 {
00477 struct tm itm;
00478 time_t ttm;
00479 const time_t s_epoch_2000 = 946684800;
00480
00481 memset(&itm, 0, sizeof(struct tm));
00482 itm.tm_year = year - 1900;
00483 itm.tm_mon = month - 1;
00484 itm.tm_mday = day;
00485 itm.tm_hour = hour;
00486 itm.tm_min = minute;
00487
00488 tzset( );
00489 ttm = mktime(&itm) + ForecastSeconds( ) - timezone - s_epoch_2000;
00490 return ttm;
00491 }
00492
00493 std::string GRIB_TIME::ValidTime( )
00494 {
00495 struct tm itm, *xtm;
00496 time_t ttm;
00497
00498 memset(&itm, 0, sizeof(struct tm));
00499 itm.tm_year = year - 1900;
00500 itm.tm_mon = month - 1;
00501 itm.tm_mday = day;
00502 itm.tm_hour = hour;
00503 itm.tm_min = minute;
00504
00505 tzset( );
00506 ttm = mktime(&itm) + ForecastSeconds( ) - timezone;
00507 xtm = gmtime(&ttm);
00508
00509 char temp[32];
00510
00511 strftime(temp, 32, "%Y-%m-%d %H:%M:00 UTC", xtm);
00512 std::string retval = temp;
00513 return retval;
00514 }
00515
00516 std::ostream& operator<< ( std::ostream& os, GRIB_TIME &t )
00517 {
00518 os << t.timestring;
00519 return os;
00520 }
00521
00522
00523
00524
00525
00526 GRIB_LEVEL::GRIB_LEVEL( )
00527 {
00528 type = GRIB_LEVEL_UNKNOWN;
00529 lv1 = 0.0;
00530 lv2 = 0.0;
00531 }
00532
00533 void GRIB_LEVEL::set(t_enum_GRIB_LEVELS type)
00534 {
00535 this->set(type, 0.0, 0.0);
00536 return;
00537 }
00538
00539 void GRIB_LEVEL::set(t_enum_GRIB_LEVELS type, float lev)
00540 {
00541 this->set(type, lev, 0.0);
00542 return;
00543 }
00544
00545 void GRIB_LEVEL::set(t_enum_GRIB_LEVELS type, float lev1, float lev2)
00546 {
00547 lv1 = 0.0;
00548 lv2 = 0.0;
00549
00550 this->type = type;
00551 switch (type)
00552 {
00553 case GRIB_LEVEL_SURFACE:
00554 case GRIB_LEVEL_CLOUD_BASE:
00555 case GRIB_LEVEL_CLOUD_TOP:
00556 case GRIB_LEVEL_ISOTHERM_0_DEG:
00557 case GRIB_LEVEL_ADIABATIC_CONDENSATION_LIFTED_FROM_SURFACE:
00558 case GRIB_LEVEL_MAXIMUM_WIND:
00559 case GRIB_LEVEL_TROPOPAUSE:
00560 case GRIB_LEVEL_NOMINAL_ATMOSPHERE_TOP:
00561 case GRIB_LEVEL_ENTIRE_ATMOSPHERE:
00562 case GRIB_LEVEL_ENTIRE_OCEAN:
00563 case GRIB_LEVEL_MEAN_SEA_LEVEL:
00564 break;
00565 case GRIB_LEVEL_ISOBARIC_mb:
00566 case GRIB_LEVEL_ALTITUDE_ABOVE_MSL_m:
00567 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_m:
00568 case GRIB_LEVEL_HYBRID:
00569 case GRIB_LEVEL_DEPTH_BELOW_SURFACE_cm:
00570 case GRIB_LEVEL_ISENTROPIC_K:
00571 case GRIB_LEVEL_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00572 case GRIB_LEVEL_POTENTIAL_VORTICITY_SURFACE_PV_UNITS:
00573 case GRIB_LEVEL_ISOBARIC_Pa:
00574 case GRIB_LEVEL_DEPTH_BELOW_SEA_m:
00575 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_HIGH_PRECISION_cm:
00576 case GRIB_LEVEL_ISOTHERMAL_K:
00577 case GRIB_LEVEL_SIGMA:
00578 case GRIB_LEVEL_ETA:
00579 lv1 = lev1;
00580 break;
00581 case GRIB_LEVEL_LAYER_HYBRID:
00582 case GRIB_LEVEL_LAYER_DEPTH_BELOW_SURFACE_cm:
00583 case GRIB_LEVEL_LAYER_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00584 case GRIB_LEVEL_LAYER_ISOBARIC_mb:
00585 case GRIB_LEVEL_LAYER_ALTITUDE_ABOVE_MSL_m:
00586 case GRIB_LEVEL_LAYER_HEIGHT_ABOVE_GROUND_m:
00587 case GRIB_LEVEL_LAYER_SIGMA:
00588 case GRIB_LEVEL_LAYER_ETA:
00589 case GRIB_LEVEL_LAYER_ISENTROPIC_K:
00590 case GRIB_LEVEL_LAYER_ISOBARIC_HIGH_PRECISION_mb:
00591 case GRIB_LEVEL_LAYER_SIGMA_HIGH_PRECISION:
00592 case GRIB_LEVEL_LAYER_ISOBARIC_MIXED_PRECISION_mb:
00593 lv1 = lev1;
00594 lv2 = lev2;
00595 break;
00596 default:
00597 std::cerr << "Unknown Level is used !!!!" << std::endl;
00598 }
00599 return;
00600 }
00601
00602 int GRIB_LEVEL::Encode( )
00603 {
00604 int retval = 0;
00605 switch (type)
00606 {
00607 case GRIB_LEVEL_SURFACE:
00608 case GRIB_LEVEL_CLOUD_BASE:
00609 case GRIB_LEVEL_CLOUD_TOP:
00610 case GRIB_LEVEL_ISOTHERM_0_DEG:
00611 case GRIB_LEVEL_ADIABATIC_CONDENSATION_LIFTED_FROM_SURFACE:
00612 case GRIB_LEVEL_MAXIMUM_WIND:
00613 case GRIB_LEVEL_TROPOPAUSE:
00614 case GRIB_LEVEL_NOMINAL_ATMOSPHERE_TOP:
00615 case GRIB_LEVEL_ENTIRE_ATMOSPHERE:
00616 case GRIB_LEVEL_ENTIRE_OCEAN:
00617 case GRIB_LEVEL_MEAN_SEA_LEVEL:
00618 break;
00619 case GRIB_LEVEL_ISOBARIC_mb:
00620 case GRIB_LEVEL_ALTITUDE_ABOVE_MSL_m:
00621 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_m:
00622 case GRIB_LEVEL_HYBRID:
00623 case GRIB_LEVEL_DEPTH_BELOW_SURFACE_cm:
00624 case GRIB_LEVEL_ISENTROPIC_K:
00625 case GRIB_LEVEL_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00626 case GRIB_LEVEL_POTENTIAL_VORTICITY_SURFACE_PV_UNITS:
00627 case GRIB_LEVEL_ISOBARIC_Pa:
00628 case GRIB_LEVEL_DEPTH_BELOW_SEA_m:
00629 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_HIGH_PRECISION_cm:
00630 retval = (int) lv1;
00631 break;
00632 case GRIB_LEVEL_ISOTHERMAL_K:
00633 retval = (int) (lv1*100.0);
00634 break;
00635 case GRIB_LEVEL_SIGMA:
00636 case GRIB_LEVEL_ETA:
00637 retval = (int) (lv1*10000.0);
00638 break;
00639 case GRIB_LEVEL_LAYER_HYBRID:
00640 case GRIB_LEVEL_LAYER_DEPTH_BELOW_SURFACE_cm:
00641 case GRIB_LEVEL_LAYER_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00642 retval = ((int) lv1) << 8 + ((int) lv2);
00643 break;
00644 case GRIB_LEVEL_LAYER_ISOBARIC_mb:
00645 retval = ((int) (lv1/10.0)) << 8 + ((int) (lv2/10.0));
00646 break;
00647 case GRIB_LEVEL_LAYER_ALTITUDE_ABOVE_MSL_m:
00648 case GRIB_LEVEL_LAYER_HEIGHT_ABOVE_GROUND_m:
00649 retval = ((int) (lv1/100.0)) << 8 + ((int) (lv2/100.0));
00650 break;
00651 case GRIB_LEVEL_LAYER_SIGMA:
00652 case GRIB_LEVEL_LAYER_ETA:
00653 retval = ((int) (lv1*100.0)) << 8 + ((int) (lv2*100.0));
00654 break;
00655 case GRIB_LEVEL_LAYER_ISENTROPIC_K:
00656 retval = ((int) (475-lv1)) << 8 + ((int) (475-lv2));
00657 break;
00658 case GRIB_LEVEL_LAYER_ISOBARIC_HIGH_PRECISION_mb:
00659 retval = ((int) (1100-lv1)) << 8 + ((int) (1100-lv2));
00660 break;
00661 case GRIB_LEVEL_LAYER_SIGMA_HIGH_PRECISION:
00662 retval = ((int) (1000*lv1)) << 8 + ((int) (1000*lv2));
00663 break;
00664 case GRIB_LEVEL_LAYER_ISOBARIC_MIXED_PRECISION_mb:
00665 retval = ((int) lv1) << 8 + ((int) (1100-lv2));
00666 break;
00667 default:
00668 std::cerr << "Unknown Level is used !!!!" << std::endl;
00669 }
00670 return retval;
00671 }
00672
00673 void GRIB_LEVEL::Decode(unsigned char *pds)
00674 {
00675 this->type = (t_enum_GRIB_LEVELS) PDS_KPDS6(pds);
00676
00677 int lev = PDS_KPDS7(pds);
00678 int o11 = lev / 256;
00679 int o12 = lev % 256;
00680
00681 lv1 = 0.0;
00682 lv2 = 0.0;
00683
00684 switch (type)
00685 {
00686 case GRIB_LEVEL_SURFACE:
00687 case GRIB_LEVEL_CLOUD_BASE:
00688 case GRIB_LEVEL_CLOUD_TOP:
00689 case GRIB_LEVEL_ISOTHERM_0_DEG:
00690 case GRIB_LEVEL_ADIABATIC_CONDENSATION_LIFTED_FROM_SURFACE:
00691 case GRIB_LEVEL_MAXIMUM_WIND:
00692 case GRIB_LEVEL_TROPOPAUSE:
00693 case GRIB_LEVEL_NOMINAL_ATMOSPHERE_TOP:
00694 case GRIB_LEVEL_ENTIRE_ATMOSPHERE:
00695 case GRIB_LEVEL_ENTIRE_OCEAN:
00696 case GRIB_LEVEL_MEAN_SEA_LEVEL:
00697 break;
00698 case GRIB_LEVEL_ISOBARIC_mb:
00699 case GRIB_LEVEL_ALTITUDE_ABOVE_MSL_m:
00700 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_m:
00701 case GRIB_LEVEL_HYBRID:
00702 case GRIB_LEVEL_DEPTH_BELOW_SURFACE_cm:
00703 case GRIB_LEVEL_ISENTROPIC_K:
00704 case GRIB_LEVEL_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00705 case GRIB_LEVEL_POTENTIAL_VORTICITY_SURFACE_PV_UNITS:
00706 case GRIB_LEVEL_ISOBARIC_Pa:
00707 case GRIB_LEVEL_DEPTH_BELOW_SEA_m:
00708 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_HIGH_PRECISION_cm:
00709 lv1 = (float) lev;
00710 break;
00711 case GRIB_LEVEL_ISOTHERMAL_K:
00712 lv1 = (float) (lev/100.0);
00713 break;
00714 case GRIB_LEVEL_SIGMA:
00715 case GRIB_LEVEL_ETA:
00716 lv1 = (float) (lev/10000.0);
00717 break;
00718 case GRIB_LEVEL_LAYER_HYBRID:
00719 case GRIB_LEVEL_LAYER_DEPTH_BELOW_SURFACE_cm:
00720 case GRIB_LEVEL_LAYER_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00721 lv1 = (float) o11;
00722 lv2 = (float) o12;
00723 break;
00724 case GRIB_LEVEL_LAYER_ISOBARIC_mb:
00725 lv1 = (float) o11*10.0;
00726 lv2 = (float) o12*10.0;
00727 break;
00728 case GRIB_LEVEL_LAYER_ALTITUDE_ABOVE_MSL_m:
00729 case GRIB_LEVEL_LAYER_HEIGHT_ABOVE_GROUND_m:
00730 lv1 = (float) o11*100.0;
00731 lv2 = (float) o12*100.0;
00732 break;
00733 case GRIB_LEVEL_LAYER_SIGMA:
00734 case GRIB_LEVEL_LAYER_ETA:
00735 lv1 = (float) o11/100.0;
00736 lv2 = (float) o12/100.0;
00737 break;
00738 case GRIB_LEVEL_LAYER_ISENTROPIC_K:
00739 lv1 = (float) (475-o11);
00740 lv2 = (float) (475-o12);
00741 break;
00742 case GRIB_LEVEL_LAYER_ISOBARIC_HIGH_PRECISION_mb:
00743 lv1 = (float) (1100-o11);
00744 lv2 = (float) (1100-o12);
00745 break;
00746 case GRIB_LEVEL_LAYER_SIGMA_HIGH_PRECISION:
00747 lv1 = (float) o11/1000.0;
00748 lv2 = (float) o12/1000.0;
00749 break;
00750 case GRIB_LEVEL_LAYER_ISOBARIC_MIXED_PRECISION_mb:
00751 lv1 = (float) o11;
00752 lv2 = (float) (1100-o12);
00753 break;
00754 default:
00755 std::cerr << "Unknown Level is used !!!!" << std::endl;
00756 }
00757 return;
00758 }
00759
00760 void GRIB_LEVEL::set_pressure(float plev)
00761 {
00762 type = GRIB_LEVEL_ISOBARIC_mb;
00763 lv1 = plev;
00764 lv2 = 0.0;
00765 return;
00766 }
00767
00768 void GRIB_LEVEL::set_height(float hgt)
00769 {
00770 type = GRIB_LEVEL_HEIGHT_ABOVE_GROUND_m;
00771 lv1 = hgt;
00772 lv2 = 0.0;
00773 return;
00774 }
00775
00776 void GRIB_LEVEL::set_seadepth(float dpt)
00777 {
00778 type = GRIB_LEVEL_DEPTH_BELOW_SEA_m;
00779 lv1 = dpt;
00780 lv2 = 0.0;
00781 return;
00782 }
00783
00784 void GRIB_LEVEL::set_isentropic(float temp)
00785 {
00786 type = GRIB_LEVEL_ISENTROPIC_K;
00787 lv1 = temp;
00788 lv2 = 0.0;
00789 return;
00790 }
00791
00792 void GRIB_LEVEL::set_sigma(float sigma)
00793 {
00794 type = GRIB_LEVEL_SIGMA;
00795 lv1 = sigma;
00796 lv2 = 0.0;
00797 return;
00798 }
00799
00800 void GRIB_LEVEL::set_surface( )
00801 {
00802 type = GRIB_LEVEL_SURFACE;
00803 lv1 = 0.0;
00804 lv2 = 0.0;
00805 return;
00806 }
00807
00808 std::ostream& operator<< ( std::ostream& os, GRIB_LEVEL &l )
00809 {
00810 os << l.leveldesc( );
00811 return os;
00812 }
00813
00814 std::string GRIB_LEVEL::leveldesc( )
00815 {
00816 std::stringstream s;
00817 switch (type)
00818 {
00819 case GRIB_LEVEL_SURFACE:
00820 s << "Ground or water surface";
00821 break;
00822 case GRIB_LEVEL_CLOUD_BASE:
00823 s << "Cloud base";
00824 break;
00825 case GRIB_LEVEL_CLOUD_TOP:
00826 s << "Cloud top";
00827 break;
00828 case GRIB_LEVEL_ISOTHERM_0_DEG:
00829 s << "Isotherm level at 0 Celsius";
00830 break;
00831 case GRIB_LEVEL_ADIABATIC_CONDENSATION_LIFTED_FROM_SURFACE:
00832 s << "Level of adiabatic condensation lifted from surface";
00833 break;
00834 case GRIB_LEVEL_MAXIMUM_WIND:
00835 s << "Level of maximum wind";
00836 break;
00837 case GRIB_LEVEL_TROPOPAUSE:
00838 s << "Tropopause level";
00839 break;
00840 case GRIB_LEVEL_NOMINAL_ATMOSPHERE_TOP:
00841 s << "Nominal atmosphere top";
00842 break;
00843 case GRIB_LEVEL_MEAN_SEA_LEVEL:
00844 s << "Mean Sea Level";
00845 break;
00846 case GRIB_LEVEL_ENTIRE_ATMOSPHERE:
00847 s << "Entire atmosphere as single layer";
00848 break;
00849 case GRIB_LEVEL_ENTIRE_OCEAN:
00850 s << "Entire ocean as single layer";
00851 break;
00852 case GRIB_LEVEL_ISOBARIC_mb:
00853 s << lv1 << " mb";
00854 break;
00855 case GRIB_LEVEL_ALTITUDE_ABOVE_MSL_m:
00856 s << lv1 << " m above MSL";
00857 break;
00858 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_m:
00859 s << lv1 << " m above ground";
00860 break;
00861 case GRIB_LEVEL_HYBRID:
00862 s << lv1 << " hybrid level";
00863 break;
00864 case GRIB_LEVEL_DEPTH_BELOW_SURFACE_cm:
00865 s << lv1 << " cm below surface";
00866 break;
00867 case GRIB_LEVEL_ISENTROPIC_K:
00868 s << lv1 << " K isentropic";
00869 break;
00870 case GRIB_LEVEL_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00871 s << lv1 << " mb difference from ground";
00872 break;
00873 case GRIB_LEVEL_POTENTIAL_VORTICITY_SURFACE_PV_UNITS:
00874 s << lv1 << " PV units";
00875 break;
00876 case GRIB_LEVEL_ISOBARIC_Pa:
00877 s << lv1 << " Pa";
00878 break;
00879 case GRIB_LEVEL_DEPTH_BELOW_SEA_m:
00880 s << lv1 << " m below sea";
00881 break;
00882 case GRIB_LEVEL_HEIGHT_ABOVE_GROUND_HIGH_PRECISION_cm:
00883 s << lv1 << " cm above ground (high precision)";
00884 break;
00885 case GRIB_LEVEL_ISOTHERMAL_K:
00886 s << lv1 << " K isotherm";
00887 break;
00888 case GRIB_LEVEL_SIGMA:
00889 s << lv1 << " sigma";
00890 break;
00891 case GRIB_LEVEL_ETA:
00892 s << lv1 << " eta";
00893 break;
00894 case GRIB_LEVEL_LAYER_ISOBARIC_mb:
00895 s << lv1 << "-" << lv2 << " mb layer";
00896 break;
00897 case GRIB_LEVEL_LAYER_ALTITUDE_ABOVE_MSL_m:
00898 s << lv1 << "-" << lv2 << " m above sea level layer";
00899 break;
00900 case GRIB_LEVEL_LAYER_HEIGHT_ABOVE_GROUND_m:
00901 s << lv1 << "-" << lv2 << " m height above ground layer";
00902 break;
00903 case GRIB_LEVEL_LAYER_SIGMA:
00904 s << lv1 << "-" << lv2 << " sigma layer";
00905 break;
00906 case GRIB_LEVEL_LAYER_ETA:
00907 s << lv1 << "-" << lv2 << " eta layer";
00908 break;
00909 case GRIB_LEVEL_LAYER_ISENTROPIC_K:
00910 s << lv1 << "-" << lv2 << " K isentropic layer";
00911 break;
00912 case GRIB_LEVEL_LAYER_HYBRID:
00913 s << lv1 << "-" << lv2 << " hybrid layer";
00914 break;
00915 case GRIB_LEVEL_LAYER_DEPTH_BELOW_SURFACE_cm:
00916 s << lv1 << "-" << lv2 << " cm below suface layer";
00917 break;
00918 case GRIB_LEVEL_LAYER_PRESSURE_DIFFERENCE_FROM_GROUND_mb:
00919 s << lv1 << "-" << lv2 << " mb difference from ground layer";
00920 break;
00921 case GRIB_LEVEL_LAYER_ISOBARIC_HIGH_PRECISION_mb:
00922 s << lv1 << "-" << lv2 << " mb high precision layer";
00923 break;
00924 case GRIB_LEVEL_LAYER_SIGMA_HIGH_PRECISION:
00925 s << lv1 << "-" << lv2 << " sigma high precision layer";
00926 break;
00927 case GRIB_LEVEL_LAYER_ISOBARIC_MIXED_PRECISION_mb:
00928 s << lv1 << "-" << lv2 << " mb mixed precision layer";
00929 break;
00930 default:
00931 throw;
00932 }
00933 return s.str( );
00934 }
00935
00936 bool GRIB_LEVEL::is_LevelPressure( )
00937 {
00938 return (type == GRIB_LEVEL_ISOBARIC_mb);
00939 }
00940
00941 bool GRIB_LEVEL::is_LevelHeight( )
00942 {
00943 return (type == GRIB_LEVEL_HEIGHT_ABOVE_GROUND_m);
00944 }
00945
00946 bool GRIB_LEVEL::is_Surface( )
00947 {
00948 return (type == GRIB_LEVEL_SURFACE);
00949 }
00950
00951 bool GRIB_LEVEL::is_MeanSeaLevel( )
00952 {
00953 return (type == GRIB_LEVEL_MEAN_SEA_LEVEL);
00954 }
00955
00956
00957
00958
00959
00960 void GRIB_GRID_regular_latlon::set(float lat1, float lat2,
00961 float lon1, float lon2,
00962 float dlat, float dlon)
00963 {
00964 this->lat1 = lat1;
00965 this->lon1 = lon1;
00966 this->lat2 = lat2;
00967 this->lon2 = lon2;
00968 this->dlat = dlat;
00969 this->dlon = dlon;
00970 return;
00971 }
00972
00973 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_regular_latlon &g )
00974 {
00975 os << " Regular latitude/longitude projection" << std::endl
00976 << " "
00977 << "Latitude " << g.lat1 << " to " << g.lat2
00978 << " by " << g.dlat << std::endl
00979 << " "
00980 << "Longitude " << g.lon1 << " to " << g.lon2
00981 << " by " << g.dlon;
00982 return os;
00983 }
00984
00985 void GRIB_GRID_gaussian_latlon::set(float lat1, float lat2,
00986 float lon1, float lon2,
00987 float dlon, float N)
00988 {
00989 this->lat1 = lat1;
00990 this->lon1 = lon1;
00991 this->lat2 = lat2;
00992 this->lon2 = lon2;
00993 this->dlon = dlon;
00994 this->N = N;
00995 return;
00996 }
00997
00998 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_gaussian_latlon &g )
00999 {
01000 os << " Gaussian latitude/longitude projection" << std::endl
01001 << " Latitude " << g.lat1 << " to " << g.lat2
01002 << " " << g.N << " latitude circles" << std::endl
01003 << " Longitude " << g.lon1 << " to " << g.lon2
01004 << " by " << g.dlon;
01005 return os;
01006 }
01007
01008 void GRIB_GRID_polar_stereographic::set(float lat1, float lon1, float lov,
01009 float dx, float dy, bool pole)
01010 {
01011 this->lat1 = lat1;
01012 this->lon1 = lon1;
01013 this->lov = lov;
01014 this->dx = dx;
01015 this->dy = dy;
01016 this->pole = pole;
01017 return;
01018 }
01019
01020 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_polar_stereographic &g )
01021 {
01022 os << " Polar stereographic projection" << std::endl
01023 << " "
01024 << "First point Latitude " << g.lat1 << ", Longitude " << g.lon1
01025 << std::endl << "Orientation longitude " << g.lov << std::endl
01026 << " "
01027 << "Dx " << g.dx << " Km, Dy " << g.dy << "Km, "
01028 << (g.pole ? "North Polar": "South Polar") << std::endl;
01029 return os;
01030 }
01031
01032 void GRIB_GRID_lambert_conformal::set(float lat1, float lon1, float lov,
01033 float latin1, float latin2,
01034 float latsp, float lonsp,
01035 float dx, float dy,
01036 bool pole, bool bipolar)
01037 {
01038 this->lat1 = lat1;
01039 this->lon1 = lon1;
01040 this->lov = lov;
01041 this->latin1 = latin1;
01042 this->latin2 = latin2;
01043 this->latsp = latsp;
01044 this->lonsp = lonsp;
01045 this->dx = dx;
01046 this->dy = dy;
01047 this->pole = pole;
01048 this->bipolar = bipolar;
01049 return;
01050 }
01051
01052 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_lambert_conformal &g )
01053 {
01054 os << " Lambert conformal projection" << std::endl
01055 << " "
01056 << "First point Latitude " << g.lat1 << ", Longitude " << g.lon1
01057 << std::endl << "Orientation longitude " << g.lov << std::endl
01058 << " "
01059 << "Secant latitude1 " << g.latin1 << ", latitude2 " << g.latin2
01060 << std::endl << "South Pole Latitude " << g.latsp
01061 << ", Longitude " << g.lonsp << std::endl
01062 << " "
01063 << "Dx " << g.dx << " Km, Dy " << g.dy << "Km, "
01064 << (g.pole ? "North Polar": "South Polar") << ", "
01065 << (g.bipolar ? " bipolar symmetric" : "one center") << std::endl;
01066 return os;
01067 }
01068
01069 void GRIB_GRID_mercator::set(float lat1, float lon1,
01070 float lat2, float lon2,
01071 float dx, float dy, float latin)
01072 {
01073 this->lat1 = lat1;
01074 this->lon1 = lon1;
01075 this->lat2 = lat2;
01076 this->lon2 = lon2;
01077 this->dx = dx;
01078 this->dy = dy;
01079 this->latin = latin;
01080 return;
01081 }
01082
01083 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_mercator &g )
01084 {
01085 os << " Mercator projection" << std::endl
01086 << " "
01087 << "Latitude " << g.lat1 << " to " << g.lat2
01088 << " by " << g.dy << " Km " << std::endl
01089 << "Longitude " << g.lon1 << " to " << g.lon2
01090 << " "
01091 << " by " << g.dx << " Km " << std::endl
01092 << " "
01093 << "Latin " << g.latin << std::endl;
01094 return os;
01095 }
01096
01097 void GRIB_GRID_spaceview::set(float lap, float lop, float dx, float dy,
01098 float Xp, float Yp, float orient, float Nr)
01099 {
01100 this->lap = lap;
01101 this->lop = lop;
01102 this->dx = dx;
01103 this->dy = dy;
01104 this->Xp = Xp;
01105 this->Yp = Yp;
01106 this->orient = orient;
01107 this->Nr = Nr;
01108 return;
01109 }
01110
01111 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_spaceview &g )
01112 {
01113 os << " Space view projection" << std::endl
01114 << " Sub satellite point lat/lon: "
01115 << g.lap << ", " << g.lop << std::endl
01116 << " Apparent diameter of earth X,Y: "
01117 << g.dx << ", " << g.dy << std::endl
01118 << " Sub satellite point X,Y: "
01119 << g.Xp << ", " << g.Yp << std::endl
01120 << " Orient " << g.orient << ", Camera H " << g.Nr
01121 << std::endl;
01122 return os;
01123 }
01124
01125 void GRIB_GRID_spherical_harmonic::set(int J, int K, int M, int rt, int csm)
01126 {
01127 this->J = J;
01128 this->K = K;
01129 this->M = M;
01130 representation = (GRIB_GRID_SPECTRAL_REPRESENTATION) rt;
01131 storage_method = (GRIB_GRID_COEFFICIENT_STORAGE_METHOD) csm;
01132 return;
01133 }
01134
01135 std::ostream& operator<< ( std::ostream& os, GRIB_GRID_spherical_harmonic &g )
01136 {
01137 os << " Spherical harmonic coefficients" << std::endl
01138 << " J, K, M : "
01139 << g.J << ", " << g.K << ", " << g.M << std::endl
01140 << " Spectral representation: "
01141 << ((g.representation == g.SR_associated_legendre_polynomials_first_kind) ?
01142 "Associated Legendre polynomials of first kind" :
01143 "Sperical armonics - complex packing")
01144 << " Storage Method : "
01145 << ((g.storage_method == g.CSM_normal) ? "Normal" : "Complex packing")
01146 << std::endl;
01147 return os;
01148 }
01149
01150 GRIB_GRID::GRIB_GRID( )
01151 {
01152 grid_code = 0;
01153 type = GRIB_GRID_UNKNOWN;
01154 is_dirincgiven = true;
01155 is_earth_spheroid = false;
01156 is_uv_grid = false;
01157 is_x_negative = false;
01158 is_y_negative = false;
01159 is_fortran = false;
01160 nxny = nx = ny = 0;
01161 }
01162
01163 void GRIB_GRID::set_size(int nx, int ny)
01164 {
01165 this->nx = nx;
01166 this->ny = ny;
01167 nxny = nx*ny;
01168 return;
01169 }
01170
01171 void GRIB_GRID::set_number(int code)
01172 {
01173 grid_code = code;
01174 return;
01175 }
01176
01177 void GRIB_GRID::set_regular_latlon(float lat1, float lat2,
01178 float lon1, float lon2,
01179 float dlat, float dlon)
01180 {
01181 grid_code = 255;
01182 type = GRIB_GRID_REGULAR_LATLON;
01183 ll.set(lat1, lat2, lon1, lon2, dlat, dlon);
01184 return;
01185 }
01186
01187 void GRIB_GRID::set_gaussian_latlon(float lat1, float lat2,
01188 float lon1, float lon2,
01189 float dlon, float N)
01190 {
01191 grid_code = 255;
01192 type = GRIB_GRID_GAUSSIAN;
01193 gl.set(lat1, lat2, lon1, lon2, dlon, N);
01194 return;
01195 }
01196
01197 void GRIB_GRID::set_polar_stereo(float lat1, float lon1, float lov,
01198 float dx, float dy, bool pole)
01199 {
01200 grid_code = 255;
01201 type = GRIB_GRID_POLAR_STEREOGRAPHIC;
01202 ps.set(lat1, lon1, lov, dx, dy, pole);
01203 return;
01204 }
01205
01206 void GRIB_GRID::set_lambert_conformal(float lat1, float lon1, float lov,
01207 float latin1, float latin2,
01208 float latsp, float lonsp,
01209 float dx, float dy,
01210 bool pole, bool bipolar)
01211 {
01212 grid_code = 255;
01213 type = GRIB_GRID_LAMBERT_CONFORMAL;
01214 lc.set(lat1, lon1, lov, latin1, latin2, latsp, lonsp, dx, dy, pole, bipolar);
01215 return;
01216 }
01217
01218 void GRIB_GRID::set_mercator(float lat1, float lon1, float lat2, float lon2,
01219 float dx, float dy, float latin)
01220 {
01221 grid_code = 255;
01222 type = GRIB_GRID_MERCATOR;
01223 mc.set(lat1, lon1, lat2, lon2, dx, dy, latin);
01224 return;
01225 }
01226
01227 void GRIB_GRID::set_spaceview(float lap, float lop, float dx, float dy,
01228 float Xp, float Yp, float orient, float Nr)
01229 {
01230 grid_code = 255;
01231 type = GRIB_GRID_SPACEVIEW;
01232 sp.set(lap, lop, dx, dy, Xp, Yp, orient, Nr);
01233 return;
01234 }
01235
01236 void GRIB_GRID::set_spherical_harmonic(int J, int K, int M, int rt, int csm)
01237 {
01238 grid_code = 255;
01239 type = GRIB_GRID_SPHERICAL_HARMONIC_COE;
01240 sa.set(J, K, M, rt, csm);
01241 return;
01242 }
01243
01244 void GRIB_GRID::set_earth_spheroid( ) { is_earth_spheroid = true; }
01245 void GRIB_GRID::set_uv_grid( ) { is_uv_grid = true; }
01246 void GRIB_GRID::set_x_negative( ) { is_x_negative = true; }
01247 void GRIB_GRID::set_y_negative( ) { is_y_negative = true; }
01248 void GRIB_GRID::set_fortran_indexing( ) { is_fortran = true; }
01249
01250 std::ostream& operator<< ( std::ostream& os, GRIB_GRID &g )
01251 {
01252 os << "Size: " << g.nx << " x " << g.ny << " ("
01253 << g.nxny << ")" << std::endl;
01254 if (g.grid_code != 0 && g.grid_code != 255)
01255 {
01256 os << " GRID IDENTIFICATION NUMBER " << g.grid_code;
01257 return os;
01258 }
01259 if (! g.is_dirincgiven)
01260 os << "Direction increment are NOT given!" << std::endl;
01261 if (g.is_earth_spheroid)
01262 os << "Earth is spheroid!" << std::endl;
01263 if (g.is_uv_grid)
01264 os << "U-V vector components are GRID relative!" << std::endl;
01265 if (g.is_x_negative)
01266 os << "X scan mode is reversed (EW)!" << std::endl;
01267 if (g.is_y_negative)
01268 os << "Y scan mode is reversed (NS)!" << std::endl;
01269 if (g.is_fortran)
01270 os << "Fortran indexing (columnar)!" << std::endl;
01271 switch (g.type)
01272 {
01273 case GRIB_GRID_REGULAR_LATLON:
01274 os << g.ll;
01275 break;
01276 case GRIB_GRID_GAUSSIAN:
01277 os << g.gl;
01278 break;
01279 case GRIB_GRID_POLAR_STEREOGRAPHIC:
01280 os << g.ps;
01281 break;
01282 case GRIB_GRID_LAMBERT_CONFORMAL:
01283 os << g.lc;
01284 break;
01285 case GRIB_GRID_MERCATOR:
01286 os << g.mc;
01287 break;
01288 case GRIB_GRID_SPACEVIEW:
01289 os << g.sp;
01290 break;
01291 case GRIB_GRID_SPHERICAL_HARMONIC_COE:
01292 os << g.sa;
01293 break;
01294 default:
01295 os << "Unsupported Grid is being used." << std::endl
01296 << "Code is: " << g.type << std::endl
01297 << "Grid Descripion is: " << std::endl
01298 << std::setw(3) << std::setfill('0')
01299 << (int) g.gds[6] << " " << (int) g.gds[7] << " "
01300 << (int) g.gds[8] << " " << (int) g.gds[9] << " "
01301 << (int) g.gds[10] << " " << (int) g.gds[11] << " "
01302 << (int) g.gds[12] << " " << (int) g.gds[13] << " "
01303 << (int) g.gds[14] << " " << (int) g.gds[15] << " "
01304 << (int) g.gds[16] << " " << (int) g.gds[17] << " "
01305 << (int) g.gds[18] << " " << (int) g.gds[19] << " "
01306 << (int) g.gds[20] << std::endl
01307 << (int) g.gds[21] << " " << (int) g.gds[22] << " "
01308 << (int) g.gds[23] << " " << (int) g.gds[24] << " "
01309 << (int) g.gds[25] << " " << (int) g.gds[26] << " "
01310 << (int) g.gds[27] << " " << (int) g.gds[28] << " "
01311 << (int) g.gds[29] << " " << (int) g.gds[30] << " "
01312 << (int) g.gds[31] << " " << (int) g.gds[32] << " "
01313 << (int) g.gds[33] << " " << (int) g.gds[34] << " "
01314 << (int) g.gds[35] << std::endl
01315 << (int) g.gds[36] << " " << (int) g.gds[37] << " "
01316 << (int) g.gds[38] << " " << (int) g.gds[39] << " "
01317 << (int) g.gds[40] << " " << (int) g.gds[41] << " "
01318 << (int) g.gds[42] << " " << (int) g.gds[43] << std::endl;
01319 }
01320 return os;
01321 }
01322
01323 unsigned char *GRIB_GRID::Encode( )
01324 {
01325 unsigned char *xgds = 0;
01326
01327 if (grid_code != 255) return 0;
01328
01329 int mode = 0, scan = 0, pole = 0;
01330 if (is_dirincgiven) mode += 128;
01331 if (is_earth_spheroid) mode += 64;
01332 if (is_uv_grid) mode += 8;
01333 if (is_x_negative) scan += 128;
01334 if (is_y_negative) scan += 64;
01335 if (is_fortran) scan += 32;
01336
01337 switch (type)
01338 {
01339 case GRIB_GRID_REGULAR_LATLON:
01340 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 32, 0,
01341 WGRIB_ENCODE_BYTE, 4, 255,
01342 WGRIB_ENCODE_BYTE, 5, GRIB_GRID_REGULAR_LATLON,
01343 WGRIB_ENCODE_2BYTES, 6, nx,
01344 WGRIB_ENCODE_2BYTES, 8, ny,
01345 WGRIB_ENCODE_S3BYTES, 10, (int) (1000.0*(ll.lat1)),
01346 WGRIB_ENCODE_S3BYTES, 13, (int) (1000.0*(ll.lon1)),
01347 WGRIB_ENCODE_BYTE, 16, mode,
01348 WGRIB_ENCODE_S3BYTES, 17, (int) (1000.0*(ll.lat2)),
01349 WGRIB_ENCODE_S3BYTES, 20, (int) (1000.0*(ll.lon2)),
01350 WGRIB_ENCODE_S2BYTES, 23, (int) (1000.0*(ll.dlon)),
01351 WGRIB_ENCODE_S2BYTES, 25, (int) (1000.0*(ll.dlat)),
01352 WGRIB_ENCODE_BYTE, 27, scan, WGRIB_ENCODE_END);
01353 break;
01354 case GRIB_GRID_POLAR_STEREOGRAPHIC:
01355 pole = ps.pole ? 0 : 128;
01356 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 32, 0,
01357 WGRIB_ENCODE_BYTE, 4, 255,
01358 WGRIB_ENCODE_BYTE, 5, GRIB_GRID_POLAR_STEREOGRAPHIC,
01359 WGRIB_ENCODE_2BYTES, 6, nx,
01360 WGRIB_ENCODE_2BYTES, 8, ny,
01361 WGRIB_ENCODE_S3BYTES, 10, (int) (1000.0*(ps.lat1)),
01362 WGRIB_ENCODE_S3BYTES, 13, (int) (1000.0*(ps.lon1)),
01363 WGRIB_ENCODE_BYTE, 16, mode,
01364 WGRIB_ENCODE_S3BYTES, 17, (int) (1000.0*(ps.lov)),
01365 WGRIB_ENCODE_S3BYTES, 20, (int) (1000.0*(ps.dx)),
01366 WGRIB_ENCODE_S3BYTES, 23, (int) (1000.0*(ps.dy)),
01367 WGRIB_ENCODE_BYTE, 26, pole,
01368 WGRIB_ENCODE_BYTE, 27, scan, WGRIB_ENCODE_END);
01369 break;
01370 case GRIB_GRID_LAMBERT_CONFORMAL:
01371 pole = 0;
01372 if ( lc.pole && lc.bipolar) pole = 64;
01373 else if (! lc.pole && ! lc.bipolar) pole = 128;
01374 else if (! lc.pole && lc.bipolar) pole = 192;
01375 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 42, 0,
01376 WGRIB_ENCODE_BYTE, 4, 255,
01377 WGRIB_ENCODE_BYTE, 5, GRIB_GRID_LAMBERT_CONFORMAL,
01378 WGRIB_ENCODE_2BYTES, 6, nx,
01379 WGRIB_ENCODE_2BYTES, 8, ny,
01380 WGRIB_ENCODE_S3BYTES, 10, (int) (1000.0*(lc.lat1)),
01381 WGRIB_ENCODE_S3BYTES, 13, (int) (1000.0*(lc.lon1)),
01382 WGRIB_ENCODE_BYTE, 16, mode,
01383 WGRIB_ENCODE_S3BYTES, 17, (int) (1000.0*(lc.lov)),
01384 WGRIB_ENCODE_S3BYTES, 20, (int) (1000.0*(lc.dx)),
01385 WGRIB_ENCODE_S3BYTES, 23, (int) (1000.0*(lc.dy)),
01386 WGRIB_ENCODE_BYTE, 26, pole,
01387 WGRIB_ENCODE_BYTE, 27, scan,
01388 WGRIB_ENCODE_S3BYTES, 28, (int) (1000.0*(lc.latin1)),
01389 WGRIB_ENCODE_S3BYTES, 31, (int) (1000.0*(lc.latin2)),
01390 WGRIB_ENCODE_S3BYTES, 34, (int) (1000.0*(lc.latsp)),
01391 WGRIB_ENCODE_S3BYTES, 37, (int) (1000.0*(lc.lonsp)),
01392 WGRIB_ENCODE_END);
01393 break;
01394 case GRIB_GRID_MERCATOR:
01395 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 42, 0,
01396 WGRIB_ENCODE_BYTE, 4, 255,
01397 WGRIB_ENCODE_BYTE, 5, GRIB_GRID_MERCATOR,
01398 WGRIB_ENCODE_2BYTES, 6, nx,
01399 WGRIB_ENCODE_2BYTES, 8, ny,
01400 WGRIB_ENCODE_S3BYTES, 10, (int) (1000.0*(mc.lat1)),
01401 WGRIB_ENCODE_S3BYTES, 13, (int) (1000.0*(mc.lon1)),
01402 WGRIB_ENCODE_BYTE, 16, mode,
01403 WGRIB_ENCODE_S3BYTES, 17, (int) (1000.0*(mc.lat2)),
01404 WGRIB_ENCODE_S3BYTES, 20, (int) (1000.0*(mc.lon2)),
01405 WGRIB_ENCODE_S3BYTES, 23, (int) (1000.0*(mc.latin)),
01406 WGRIB_ENCODE_BYTE, 27, scan,
01407 WGRIB_ENCODE_S3BYTES, 28, (int) (1000.0*(mc.dx)),
01408 WGRIB_ENCODE_S3BYTES, 31, (int) (1000.0*(mc.dy)),
01409 WGRIB_ENCODE_END);
01410 break;
01411 case GRIB_GRID_SPACEVIEW:
01412 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 44, 0,
01413 WGRIB_ENCODE_BYTE, 4, 255,
01414 WGRIB_ENCODE_BYTE, 5, GRIB_GRID_SPACEVIEW,
01415 WGRIB_ENCODE_2BYTES, 6, nx,
01416 WGRIB_ENCODE_2BYTES, 8, ny,
01417 WGRIB_ENCODE_S3BYTES, 10, (int) (1000.0*(sp.lap)),
01418 WGRIB_ENCODE_S3BYTES, 13, (int) (1000.0*(sp.lop)),
01419 WGRIB_ENCODE_BYTE, 16, mode,
01420 WGRIB_ENCODE_S3BYTES, 17, (int) (1000.0*(sp.dx)),
01421 WGRIB_ENCODE_S3BYTES, 20, (int) (1000.0*(sp.dy)),
01422 WGRIB_ENCODE_S2BYTES, 23, (int) sp.Xp,
01423 WGRIB_ENCODE_S2BYTES, 25, (int) sp.Yp,
01424 WGRIB_ENCODE_BYTE, 27, scan,
01425 WGRIB_ENCODE_S3BYTES, 28, (int) (1000.0*(sp.orient)),
01426 WGRIB_ENCODE_S3BYTES, 31, (int) (1000.0*(sp.Nr)),
01427 WGRIB_ENCODE_END);
01428 break;
01429 case GRIB_GRID_SPHERICAL_HARMONIC_COE:
01430 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 32, 0,
01431 WGRIB_ENCODE_BYTE, 4, 255,
01432 WGRIB_ENCODE_BYTE, 5, GRIB_GRID_SPHERICAL_HARMONIC_COE,
01433 WGRIB_ENCODE_2BYTES, 6, sa.J,
01434 WGRIB_ENCODE_2BYTES, 8, sa.K,
01435 WGRIB_ENCODE_2BYTES, 10, sa.M,
01436 WGRIB_ENCODE_BYTE, 12, sa.representation,
01437 WGRIB_ENCODE_BYTE, 13, sa.storage_method,
01438 WGRIB_ENCODE_END);
01439 break;
01440 default:
01441 xgds = mk_GDS(NULL, WGRIB_ENCODE_INIT, 32, 0,
01442 WGRIB_ENCODE_BYTE, 4, 255,
01443 WGRIB_ENCODE_END);
01444 memcpy(xgds+4, gds+4, 28);
01445 break;
01446 }
01447 return xgds;
01448 }
01449
01450 void GRIB_GRID::Decode(unsigned char *pds,
01451 unsigned char *xgds,
01452 unsigned char *bms,
01453 unsigned char *bds)
01454 {
01455
01456 ny = 1;
01457
01458 if (xgds != NULL)
01459 {
01460 if (GDS_LEN(xgds) > 44)
01461 {
01462 std::cerr << "Sorry, thinned grids NOT supported." << std::endl;
01463 throw;
01464 }
01465
01466 memcpy(gds, xgds, GDS_LEN(xgds));
01467
01468 GDS_grid(gds, bds, &nx, &ny, &nxny);
01469 type = (t_enum_GRIB_GRIDS) GDS_DataType(gds);
01470 switch (type)
01471 {
01472 case GRIB_GRID_REGULAR_LATLON:
01473 ll.set(0.001*GDS_LatLon_La1(gds), 0.001*GDS_LatLon_La2(gds),
01474 0.001*GDS_LatLon_Lo1(gds), 0.001*GDS_LatLon_Lo2(gds),
01475 0.001*GDS_LatLon_dy(gds), 0.001*GDS_LatLon_dx(gds));
01476 break;
01477 case GRIB_GRID_GAUSSIAN:
01478 gl.set(0.001*GDS_LatLon_La1(gds), 0.001*GDS_LatLon_La2(gds),
01479 0.001*GDS_LatLon_Lo1(gds), 0.001*GDS_LatLon_Lo2(gds),
01480 0.001*GDS_LatLon_dx(gds), 0.001*GDS_Gaussian_nlat(gds));
01481 break;
01482 case GRIB_GRID_POLAR_STEREOGRAPHIC:
01483 ps.set(0.001*GDS_Polar_La1(gds), 0.001*GDS_Polar_Lo1(gds),
01484 0.001*GDS_Polar_Lov(gds),
01485 0.001*GDS_Polar_Dx(gds), 0.001*GDS_Polar_Dy(gds),
01486 (GDS_Polar_pole(gds) ? true : false));
01487 break;
01488 case GRIB_GRID_LAMBERT_CONFORMAL:
01489 lc.set(0.001*GDS_Lambert_La1(gds), 0.001*GDS_Lambert_Lo1(gds),
01490 0.001*GDS_Lambert_Lov(gds), 0.001*GDS_Lambert_Latin1(gds),
01491 0.001*GDS_Lambert_Latin2(gds), 0.001*GDS_Lambert_LatSP(gds),
01492 0.001*GDS_Lambert_LonSP(gds), 0.001*GDS_Lambert_dx(gds),
01493 0.001*GDS_Lambert_dy(gds), (GDS_Lambert_NP(gds) ? true : false),
01494 (((gds[26] & 64) == 0) ? true : false));
01495 break;
01496 case GRIB_GRID_MERCATOR:
01497 mc.set(0.001*GDS_Merc_La1(gds), 0.001*GDS_Merc_Lo1(gds),
01498 0.001*GDS_Merc_La2(gds), 0.001*GDS_Merc_Lo2(gds),
01499 0.001*GDS_Merc_dx(gds), 0.001*GDS_Merc_dy(gds),
01500 0.001*GDS_Merc_Latin(gds));
01501 break;
01502 case GRIB_GRID_SPHERICAL_HARMONIC_COE:
01503 sa.set(GDS_Harmonic_nj(gds), GDS_Harmonic_nk(gds),
01504 GDS_Harmonic_nm(gds), GDS_Harmonic_type(gds),
01505 GDS_Harmonic_mode(gds));
01506 break;
01507 case GRIB_GRID_SPACEVIEW:
01508 sp.set(0.001*GDS_Spaceview_lap(gds), 0.001*GDS_Spaceview_lop(gds),
01509 0.001*GDS_Spaceview_dx(gds), 0.001*GDS_Spaceview_dy(gds),
01510 GDS_Spaceview_Xp(gds), GDS_Spaceview_Yp(gds),
01511 0.001*GDS_Spaceview_or(gds), 0.001*GDS_Spaceview_nr(gds));
01512 break;
01513 default:
01514 break;
01515 }
01516
01517 is_dirincgiven = (gds[16] & 128);
01518 is_earth_spheroid = (gds[16] & 64);
01519 is_uv_grid = (gds[16] & 8);
01520 is_x_negative = (gds[27] & 128);
01521 is_y_negative = (gds[27] & 64);
01522 is_fortran = (gds[27] & 32);
01523 }
01524 else if (bms != NULL)
01525 nxny = nx = ((bms[4]<<8) + bms[5]);
01526 else
01527 {
01528 if (BDS_NumBits(bds) == 0)
01529 nxny = nx = 1;
01530 else
01531 nxny = nx = BDS_NValues(bds);
01532 }
01533
01534 grid_code = PDS_Grid(pds);
01535
01536 return;
01537 }
01538
01539
01540
01541
01542
01543 GRIB_FIELD::GRIB_FIELD( )
01544 {
01545 varcode = GRIB_PARAMETER_UNKNOWN;
01546 center = GRIB_CENTER_LOCAL;
01547 subcenter = GRIB_SUBCENTER_LOCAL;
01548 table = GRIB_TABLE_INTERNATIONAL;
01549 process = GRIB_PROCESS_LOCAL;
01550 vals = 0;
01551 size = 0;
01552 varname = "Unknown";
01553 varunit = "Unknown";
01554 vardesc = "Unknown";
01555 numbits = 0;
01556 refvalue = 0.0;
01557 decimalscale = 0;
01558 binscale = 0;
01559 undef_high = 9.9991e20;
01560 undef_low = 9.9991e20;
01561 }
01562
01563 GRIB_FIELD::~GRIB_FIELD( )
01564 {
01565 if (vals) delete [ ] vals;
01566 }
01567
01568 void GRIB_FIELD::Decode(unsigned char *pds, unsigned char *bds,
01569 unsigned char *bms, size_t nxny)
01570 {
01571 if (size != nxny)
01572 {
01573 if (vals) delete [ ] vals;
01574 vals = new float[nxny];
01575 }
01576
01577 if (vals == 0)
01578 vals = new float[nxny];
01579
01580 size = nxny;
01581
01582 center = PDS_Center(pds);
01583 subcenter = PDS_Subcenter(pds);
01584 table = PDS_Vsn(pds);
01585 process = PDS_Model(pds);
01586
01587 varcode = PDS_PARAM(pds);
01588
01589 varname = (Parm_Table(center, subcenter, table, process) + varcode)->name;
01590 vardesc = (Parm_Table(center, subcenter, table, process) + varcode)->comment;
01591 size_t indx0 = vardesc.find_first_of('[', 0);
01592 size_t indx1 = vardesc.find_first_of(']', 0);
01593 varunit = vardesc.substr(indx0+1,indx1-indx0-1);
01594
01595 numbits = BDS_NumBits(bds);
01596 decimalscale = PDS_DecimalScale(pds);
01597 refvalue = BDS_RefValue(bds);
01598 binscale = BDS_BinScale(bds);
01599
01600 double temp = int_power(10.0, - decimalscale);
01601
01602 BDS_unpack(vals, bds, BMS_bitmap(bms), numbits, size,
01603 temp * refvalue, temp * int_power(2.0, binscale));
01604
01605 return;
01606 }
01607
01608 std::string GRIB_FIELD::VarName( ) { return varname; }
01609 std::string GRIB_FIELD::VarUnit( ) { return varunit; }
01610 std::string GRIB_FIELD::VarDescription( ) { return vardesc; }
01611
01612 void GRIB_FIELD::set_table(int c, int s, int t, int p)
01613 {
01614 center = c;
01615 subcenter = s;
01616 table = t;
01617 process = p;
01618 return;
01619 }
01620
01621 void GRIB_FIELD::set_scale(int scale)
01622 {
01623 decimalscale = scale;
01624 return;
01625 }
01626
01627 void GRIB_FIELD::set_field(int v, float *values, size_t s, float uh, float ul)
01628 {
01629 varcode = v;
01630 varname = (Parm_Table(center, subcenter, table, process) + varcode)->name;
01631 vardesc = (Parm_Table(center, subcenter, table, process) + varcode)->comment;
01632 size_t indx0 = vardesc.find_first_of('[', 0);
01633 size_t indx1 = vardesc.find_first_of(']', 0);
01634 varunit = vardesc.substr(indx0+1,indx1-indx0-1);
01635 if (size != s || vals == 0)
01636 {
01637 size = s;
01638 if (vals) delete [ ] vals;
01639 vals = new float[s];
01640 }
01641 memcpy(vals, values, s*sizeof(float));
01642 undef_high = uh;
01643 undef_low = ul;
01644 return;
01645 }
01646
01647
01648
01649
01650
01651 GRIB_MESSAGE::GRIB_MESSAGE( )
01652 {
01653 message = 0;
01654 pds = 0;
01655 gds = 0;
01656 bms = 0;
01657 bds = 0;
01658 reclen = 0;
01659 }
01660
01661 GRIB_MESSAGE::~GRIB_MESSAGE( ) { if (message) delete [ ] message; }
01662
01663 int GRIB_MESSAGE::Decode(unsigned char *msg, size_t msgsize)
01664 {
01665 if (message) delete [ ] message;
01666
01667 unsigned char *pnt;
01668
01669 message = new unsigned char[msgsize];
01670 memcpy(message, msg, msgsize);
01671 reclen = msgsize;
01672
01673 pds = message + 8;
01674 pnt = pds + PDS_LEN(pds);
01675 if (PDS_HAS_GDS(pds))
01676 {
01677 gds = pnt;
01678 pnt += GDS_LEN(gds);
01679 }
01680 if (PDS_HAS_BMS(pds))
01681 {
01682 bms = pnt;
01683 pnt += BMS_LEN(bms);
01684 }
01685 bds = pnt;
01686 pnt += BDS_LEN(bds);
01687
01688 if ((size_t) (pnt-message+4) != reclen)
01689 {
01690 std::cerr << "Message GRIB lenght is inconsistent" << std::endl;
01691 return -1;
01692 }
01693
01694 gtime.Decode(pds);
01695 level.Decode(pds);
01696 grid.Decode(pds, gds, bms, bds);
01697 field.Decode(pds, bds, bms, (size_t) grid.nxny);
01698
01699 return 0;
01700 }
01701
01702 void GRIB_MESSAGE::Encode( )
01703 {
01704 static unsigned char header[8] = {'G', 'R', 'I', 'B', ' ', ' ', ' ', '\1'};
01705 static unsigned char trailer[4] = {'7', '7', '7', '7'};
01706
01707 unsigned char *tenc = gtime.Encode( );
01708 size_t pdslen = 0, gdslen = 0, bmslen = 0, bdslen = 0;
01709
01710 if (grid.grid_code == 0)
01711 {
01712 std::cerr << "Undefined GRID !" << std::endl;
01713 throw;
01714 }
01715
01716 if (level.type == GRIB_LEVEL_UNKNOWN)
01717 {
01718 std::cerr << "Undefined LEVEL !" << std::endl;
01719 throw;
01720 }
01721
01722 if (gtime.timeunit == GRIB_TIMEUNIT_UNKNOWN)
01723 {
01724 std::cerr << "Undefined TIME !" << std::endl;
01725 throw;
01726 }
01727
01728 if (field.varcode == GRIB_PARAMETER_UNKNOWN)
01729 {
01730 std::cerr << "Undefined FIELD !" << std::endl;
01731 throw;
01732 }
01733
01734 pds = mk_PDS(NULL, WGRIB_ENCODE_INIT, 28, 0,
01735 WGRIB_ENCODE_BYTE, 3, field.table,
01736 WGRIB_ENCODE_BYTE, 4, field.center,
01737 WGRIB_ENCODE_BYTE, 5, field.process,
01738 WGRIB_ENCODE_BYTE, 6, grid.grid_code,
01739 WGRIB_ENCODE_BYTE, 8, field.varcode,
01740 WGRIB_ENCODE_BYTE, 9, level.type,
01741 WGRIB_ENCODE_2BYTES, 10, level.Encode( ),
01742 WGRIB_ENCODE_BYTE, 12, tenc[0],
01743 WGRIB_ENCODE_BYTE, 13, tenc[1],
01744 WGRIB_ENCODE_BYTE, 14, tenc[2],
01745 WGRIB_ENCODE_BYTE, 15, tenc[3],
01746 WGRIB_ENCODE_BYTE, 16, tenc[4],
01747 WGRIB_ENCODE_BYTE, 17, tenc[5],
01748 WGRIB_ENCODE_BYTE, 18, tenc[6],
01749 WGRIB_ENCODE_BYTE, 19, tenc[7],
01750 WGRIB_ENCODE_BYTE, 20, tenc[8],
01751 WGRIB_ENCODE_BYTE, 21, tenc[9],
01752 WGRIB_ENCODE_BYTE, 22, tenc[10],
01753 WGRIB_ENCODE_BYTE, 23, tenc[11],
01754 WGRIB_ENCODE_BYTE, 24, tenc[12],
01755 WGRIB_ENCODE_BYTE, 25, field.subcenter,
01756 WGRIB_ENCODE_S2BYTES, 26, field.decimalscale,
01757 WGRIB_ENCODE_END);
01758
01759 pdslen = PDS_LEN(pds);
01760
01761 gds = grid.Encode( );
01762 if (gds)
01763 {
01764 pds[7] |= 128;
01765 gdslen = GDS_LEN(gds);
01766 }
01767
01768 int mask_nxny = field.size;
01769 float *tocomp = new float[field.size];
01770 memcpy(tocomp, field.vals, field.size*sizeof(float));
01771
01772
01773 bms = mk_BMS(pds, tocomp, &mask_nxny, field.undef_low, field.undef_high);
01774 if (bms)
01775 bmslen = BMS_LEN(bms);
01776
01777 bds = mk_BDS(pds, tocomp, mask_nxny);
01778 bdslen = BDS_LEN(bds);
01779
01780 delete [ ] tocomp;
01781
01782 reclen = pdslen + gdslen + bmslen + bdslen + 12;
01783
01784 header[4] = (reclen >> 16) & 255;
01785 header[5] = (reclen >> 8) & 255;
01786 header[6] = (reclen ) & 255;
01787
01788 if (message) delete [ ] message;
01789 message = new unsigned char[reclen];
01790
01791 size_t lpos = 0;
01792 memcpy(message, header, 8*sizeof(unsigned char));
01793 lpos += 8;
01794 memcpy(message+lpos, pds, pdslen*sizeof(unsigned char));
01795 free(pds);
01796 pds = message+lpos;
01797 lpos += pdslen;
01798 if (gds)
01799 {
01800 memcpy(message+lpos, gds, gdslen*sizeof(unsigned char));
01801 free(gds);
01802 gds = message+lpos;
01803 lpos += gdslen;
01804 }
01805 if (bms)
01806 {
01807 memcpy(message+lpos, bms, bmslen*sizeof(unsigned char));
01808 free(bms);
01809 bms = message+lpos;
01810 lpos += bmslen;
01811 }
01812 memcpy(message+lpos, bds, bdslen*sizeof(unsigned char));
01813 free(bds);
01814 bds = message+lpos;
01815 lpos += bdslen;
01816
01817 memcpy(message+lpos, trailer, 4*sizeof(unsigned char));
01818 lpos += 4;
01819
01820 if (reclen != lpos) throw;
01821
01822 return;
01823 }
01824
01825 void GRIB_MESSAGE::set_grid(GRIB_GRID &g)
01826 {
01827 grid.grid_code = g.grid_code;
01828 grid.type = g.type;
01829 grid.nxny = g.nxny;
01830 grid.nx = g.nx;
01831 grid.ny = g.ny;
01832 grid.is_dirincgiven = g.is_dirincgiven;
01833 grid.is_earth_spheroid = g.is_earth_spheroid;
01834 grid.is_uv_grid = g.is_uv_grid;
01835 grid.is_x_negative = g.is_x_negative;
01836 grid.is_y_negative = g.is_y_negative;
01837 grid.is_fortran = g.is_fortran;
01838 switch (g.type)
01839 {
01840 case GRIB_GRID_REGULAR_LATLON:
01841 grid.ll.lat1 = g.ll.lat1;
01842 grid.ll.lon1 = g.ll.lon1;
01843 grid.ll.lat2 = g.ll.lat2;
01844 grid.ll.lon2 = g.ll.lon2;
01845 grid.ll.dlat = g.ll.dlat;
01846 grid.ll.dlon = g.ll.dlon;
01847 break;
01848 case GRIB_GRID_GAUSSIAN:
01849 grid.gl.lat1 = g.gl.lat1;
01850 grid.gl.lon1 = g.gl.lon1;
01851 grid.gl.lat2 = g.gl.lat2;
01852 grid.gl.lon2 = g.gl.lon2;
01853 grid.gl.dlon = g.gl.dlon;
01854 grid.gl.N = g.gl.N;
01855 break;
01856 case GRIB_GRID_POLAR_STEREOGRAPHIC:
01857 grid.ps.lat1 = g.ps.lat1;
01858 grid.ps.lon1 = g.ps.lon1;
01859 grid.ps.lov = g.ps.lov;
01860 grid.ps.dx = g.ps.dx;
01861 grid.ps.dy = g.ps.dy;
01862 grid.ps.pole = g.ps.pole;
01863 break;
01864 case GRIB_GRID_LAMBERT_CONFORMAL:
01865 grid.lc.lat1 = g.lc.lat1;
01866 grid.lc.lon1 = g.lc.lon1;
01867 grid.lc.lov = g.lc.lov;
01868 grid.lc.latin1 = g.lc.latin1;
01869 grid.lc.latin2 = g.lc.latin2;
01870 grid.lc.latsp = g.lc.latsp;
01871 grid.lc.lonsp = g.lc.lonsp;
01872 grid.lc.dx = g.lc.dx;
01873 grid.lc.dy = g.lc.dy;
01874 grid.lc.pole = g.lc.pole;
01875 grid.lc.bipolar = g.lc.bipolar;
01876 break;
01877 case GRIB_GRID_MERCATOR:
01878 grid.mc.lat1 = g.mc.lat1;
01879 grid.mc.lon1 = g.mc.lon1;
01880 grid.mc.lat2 = g.mc.lat2;
01881 grid.mc.lon2 = g.mc.lon2;
01882 grid.mc.dx = g.mc.dx;
01883 grid.mc.dy = g.mc.dy;
01884 grid.mc.latin = g.mc.latin;
01885 break;
01886 case GRIB_GRID_SPHERICAL_HARMONIC_COE:
01887 grid.sa.J = g.sa.J;
01888 grid.sa.K = g.sa.K;
01889 grid.sa.M = g.sa.M;
01890 grid.sa.representation = g.sa.representation;
01891 grid.sa.storage_method = g.sa.storage_method;
01892 break;
01893 case GRIB_GRID_SPACEVIEW:
01894 grid.sp.lap = g.sp.lap;
01895 grid.sp.lop = g.sp.lop;
01896 grid.sp.dx = g.sp.dx;
01897 grid.sp.dy = g.sp.dy;
01898 grid.sp.Xp = g.sp.Xp;
01899 grid.sp.Yp = g.sp.Yp;
01900 grid.sp.orient = g.sp.orient;
01901 grid.sp.Nr = g.sp.Nr;
01902 break;
01903 default:
01904 break;
01905 }
01906 return;
01907 }
01908
01909 void GRIB_MESSAGE::set_time(GRIB_TIME &t)
01910 {
01911 gtime.reftime = t.reftime;
01912 gtime.timestring = t.timestring;
01913 gtime.year = t.year;
01914 gtime.month = t.month;
01915 gtime.day = t.day;
01916 gtime.hour = t.hour;
01917 gtime.minute = t.minute;
01918 gtime.p1 = t.p1;
01919 gtime.p2 = t.p2;
01920 gtime.p1p2 = t.p1p2;
01921 gtime.naveraged = t.naveraged;
01922 gtime.nmissing = t.nmissing;
01923 gtime.timeunit = t.timeunit;
01924 gtime.timerange = t.timerange;
01925 return;
01926 }
01927
01928 void GRIB_MESSAGE::set_level(GRIB_LEVEL &l)
01929 {
01930 level.type = l.type;
01931 level.lv1 = l.lv1;
01932 level.lv2 = l.lv2;
01933 return;
01934 }
01935
01936 void GRIB_MESSAGE::set_field(GRIB_FIELD &f)
01937 {
01938 if (field.vals) delete [ ] field.vals;
01939 field.vals = 0;
01940 field.size = 0;
01941
01942 if (f.size)
01943 {
01944 field.size = f.size;
01945 field.vals = new float[f.size];
01946 memcpy(field.vals, f.vals, f.size*sizeof(float));
01947 }
01948 else
01949 field.vals = 0;
01950
01951 field.varcode = f.varcode;
01952 field.varname = f.varname;
01953 field.varunit = f.varunit;
01954 field.vardesc = f.vardesc;
01955 field.numbits = f.numbits;
01956 field.refvalue = f.refvalue;
01957 field.decimalscale = f.decimalscale;
01958 field.binscale = f.binscale;
01959 field.undef_high = f.undef_high;
01960 field.undef_low = f.undef_low;
01961 return;
01962 }
01963
01964
01965
01966
01967
01968
01969 GRIB_FILE::GRIB_FILE( )
01970 {
01971 pos = 0;
01972 size = 0;
01973 limit = 0;
01974 data = 0;
01975 fd = -1;
01976 }
01977
01978 int GRIB_FILE::OpenRead(char *fname)
01979 {
01980 struct stat finfo;
01981
01982 if (stat(fname, &finfo) < 0)
01983 {
01984 std::cerr << "Error opening input file " << fname << std::endl;
01985 return -1;
01986 }
01987
01988 size = finfo.st_size;
01989 limit = size-4;
01990
01991 fd = open(fname, O_RDONLY);
01992 if (fd < 0)
01993 {
01994 std::cerr << "Error opening input file " << fname << std::endl;
01995 return -1;
01996 }
01997
01998 data = (unsigned char *) mmap(0, size, PROT_READ, MAP_PRIVATE, fd, 0);
01999 if (data == MAP_FAILED)
02000 {
02001 std::cerr << "Error reading input file " << fname << std::endl;
02002 return -1;
02003 }
02004
02005 return 0;
02006 }
02007
02008 int GRIB_FILE::OpenRead(std::string fname)
02009 {
02010 return this->OpenRead(fname.c_str( ));
02011 }
02012
02013 int GRIB_FILE::OpenWrite(char *fname)
02014 {
02015 fd = open(fname, O_WRONLY|O_CREAT|O_TRUNC, S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH);
02016 if (fd < 0)
02017 {
02018 std::cerr << "Error opening output file " << fname << std::endl;
02019 return -1;
02020 }
02021
02022 return 0;
02023 }
02024
02025 int GRIB_FILE::OpenWrite(std::string fname)
02026 {
02027 return this->OpenWrite(fname.c_str( ));
02028 }
02029
02030 int GRIB_FILE::Append(char *fname)
02031 {
02032 fd = open(fname, O_WRONLY|O_CREAT|O_APPEND, S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH);
02033 if (fd < 0)
02034 {
02035 std::cerr << "Error opening output file " << fname << std::endl;
02036 return -1;
02037 }
02038
02039 return 0;
02040 }
02041
02042 int GRIB_FILE::Append(std::string fname)
02043 {
02044 return this->Append(fname.c_str( ));
02045 }
02046
02047 int GRIB_FILE::Close( )
02048 {
02049 if (data) munmap(data, size);
02050 if (fd > 0) close(fd);
02051 pos = 0;
02052 size = 0;
02053 limit = 0;
02054 data = 0;
02055 fd = -1;
02056 return 0;
02057 }
02058
02059 GRIB_FILE::~GRIB_FILE( ) { (void) this->Close( ); }
02060
02061 int GRIB_FILE::ReadMessage(GRIB_MESSAGE &message)
02062 {
02063 size_t reclen = 0;
02064
02065 if (pos >= limit)
02066 return 1;
02067
02068 char *cpnt = (char *) data;
02069
02070 while (pos < limit)
02071 {
02072 if (*(cpnt+pos) == 'G' &&
02073 *(cpnt+pos+1) == 'R' &&
02074 *(cpnt+pos+2) == 'I' &&
02075 *(cpnt+pos+3) == 'B') break;
02076 ++pos;
02077 }
02078
02079 if (pos >= limit)
02080 return 1;
02081
02082 unsigned char *msg = data+pos;
02083
02084 reclen = (msg[4]<<16) + (msg[5]<<8) + msg[6];
02085
02086 if (msg[7] != 1)
02087 {
02088 std::cerr << "Message GRIB version is not 1: " << (int) msg[7]
02089 << std::endl;
02090 return -1;
02091 }
02092
02093 pos = pos + reclen;
02094 return message.Decode(msg, reclen);
02095 }
02096
02097 int GRIB_FILE::WriteMessage(GRIB_MESSAGE &msg)
02098 {
02099 size_t wbytes;
02100
02101 msg.Encode( );
02102 if (msg.message == 0) throw;
02103
02104 wbytes = write(fd, (void *) msg.message, msg.reclen);
02105 if (wbytes != msg.reclen)
02106 {
02107 std::cerr << "Error writing output message !" << std::endl;
02108 throw;
02109 }
02110
02111 return 0;
02112 }
02113
02114 #ifdef TESTME
02115
02116 int main(int argc, char *argv[])
02117 {
02118 if (argc < 2) return -1;
02119
02120 GRIB_FILE gf;
02121 GRIB_FILE gn;
02122
02123 int ret = gf.OpenRead(argv[1]);
02124 if (ret != 0) return -1;
02125
02126 GRIB_MESSAGE m;
02127
02128 while ((ret = gf.ReadMessage(m)) == 0)
02129 {
02130 std::cout << "Message: " << m.field.VarDescription( ) << std::endl;
02131 std::cout << "Grid : " << m.grid << std::endl
02132 << "Level: " << m.level << std::endl
02133 << "Time : " << m.gtime << std::endl << std::endl;
02134
02135 }
02136
02137 if (ret < 0) return -1;
02138
02139 ret = gf.Close( );
02140 if (ret != 0) return -1;
02141
02142 GRIB_MESSAGE m1;
02143 GRIB_GRID g;
02144 GRIB_LEVEL l;
02145 GRIB_TIME t;
02146 GRIB_FIELD f;
02147
02148 ret = gn.OpenWrite("outtest.grb");
02149 if (ret != 0) return -1;
02150
02151 float *ff = new float[65160];
02152 for (int i = 0; i < 65160; i ++)
02153 ff[i] = 273.15 + i/360;
02154
02155
02156 t.set_referencetime(2004, 7, 3, 12, 0);
02157 t.set_forecast_hour(12);
02158
02159
02160 l.set_pressure(1000.0);
02161
02162
02163 g.set_size(360, 181);
02164 g.set_regular_latlon(90.0, -90.0, 0.0, -1.0, 1.0, 1.0);
02165
02166
02167 f.set_field(GRIB_PARAMETER_TMP__, ff, 65160, 9.9991e20, 9.9991e20);
02168
02169 m1.set_time(t);
02170 m1.set_grid(g);
02171 m1.set_level(l);
02172 m1.set_field(f);
02173
02174 gn.WriteMessage(m1);
02175
02176 for (int i = 0; i < 65160; i ++)
02177 ff[i] = 10000 + i/360;
02178
02179 f.set_field(GRIB_PARAMETER_PRES_, ff, 65160, 9.9991e20, 9.9991e20);
02180
02181 m1.set_field(f);
02182
02183 gn.WriteMessage(m1);
02184
02185 ret = gn.Close( );
02186 if (ret != 0) return -1;
02187
02188 delete [ ] ff;
02189
02190 std::cout << "Done." << std::endl;
02191 return 0;
02192 }
02193
02194 #endif
