adj_cellhd.c 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894
  1. /*!
  2. * \file lib/gis/adj_cellhd.c
  3. *
  4. * \brief GIS Library - CELL header adjustment.
  5. *
  6. * (C) 2001-2009 by the GRASS Development Team
  7. *
  8. * This program is free software under the GNU General Public License
  9. * (>=v2). Read the file COPYING that comes with GRASS for details.
  10. *
  11. * \author Original author CERL
  12. */
  13. #include <math.h>
  14. #include <string.h>
  15. #include <grass/gis.h>
  16. #include <grass/glocale.h>
  17. #define LL_TOLERANCE 10
  18. /* TODO: find good thresholds */
  19. /* deviation measured in cells */
  20. static double llepsilon = 0.01;
  21. static double fpepsilon = 1.0e-9;
  22. static int ll_wrap(struct Cell_head *cellhd);
  23. static int ll_check_ns(struct Cell_head *cellhd);
  24. static int ll_check_ew(struct Cell_head *cellhd);
  25. /*!
  26. * \brief Adjust cell header.
  27. *
  28. * This function fills in missing parts of the input cell header (or
  29. * region). It also makes projection-specific adjustments. The
  30. * <i>cellhd</i> structure must have its <i>north, south, east,
  31. * west</i>, and <i>proj</i> fields set.
  32. *
  33. * If <i>row_flag</i> is true, then the north-south resolution is
  34. * computed from the number of <i>rows</i> in the <i>cellhd</i>
  35. * structure. Otherwise the number of <i>rows</i> is computed from the
  36. * north-south resolution in the structure, similarly for
  37. * <i>col_flag</i> and the number of columns and the east-west
  38. * resolution.
  39. *
  40. * <b>Note:</b> 3D values are not adjusted.
  41. *
  42. * \param[in,out] cellhd pointer to Cell_head structure
  43. * \param row_flag compute n-s resolution
  44. * \param col_flag compute e-w resolution
  45. */
  46. void G_adjust_Cell_head(struct Cell_head *cellhd, int row_flag, int col_flag)
  47. {
  48. double old_res;
  49. if (!row_flag) {
  50. if (cellhd->ns_res <= 0)
  51. G_fatal_error(_("Illegal n-s resolution value: %g"),
  52. cellhd->ns_res);
  53. }
  54. else {
  55. if (cellhd->rows <= 0)
  56. G_fatal_error(_("Illegal number of rows: %d"
  57. " (resolution is %g)"),
  58. cellhd->rows, cellhd->ns_res);
  59. }
  60. if (!col_flag) {
  61. if (cellhd->ew_res <= 0)
  62. G_fatal_error(_("Illegal e-w resolution value: %g"),
  63. cellhd->ew_res);
  64. }
  65. else {
  66. if (cellhd->cols <= 0)
  67. G_fatal_error(_("Illegal number of columns: %d"
  68. " (resolution is %g)"),
  69. cellhd->cols, cellhd->ew_res);
  70. }
  71. /* check the edge values */
  72. if (cellhd->north <= cellhd->south) {
  73. if (cellhd->proj == PROJECTION_LL)
  74. G_fatal_error(_("North must be north of South,"
  75. " but %g (north) <= %g (south"),
  76. cellhd->north, cellhd->south);
  77. else
  78. G_fatal_error(_("North must be larger than South,"
  79. " but %g (north) <= %g (south"),
  80. cellhd->north, cellhd->south);
  81. }
  82. ll_wrap(cellhd);
  83. if (cellhd->east <= cellhd->west)
  84. G_fatal_error(_("East must be larger than West,"
  85. " but %g (east) <= %g (west)"),
  86. cellhd->east, cellhd->west);
  87. /* compute rows and columns, if not set */
  88. if (!row_flag) {
  89. cellhd->rows =
  90. (cellhd->north - cellhd->south +
  91. cellhd->ns_res / 2.0) / cellhd->ns_res;
  92. if (cellhd->rows == 0)
  93. cellhd->rows = 1;
  94. }
  95. if (!col_flag) {
  96. cellhd->cols =
  97. (cellhd->east - cellhd->west +
  98. cellhd->ew_res / 2.0) / cellhd->ew_res;
  99. if (cellhd->cols == 0)
  100. cellhd->cols = 1;
  101. }
  102. if (cellhd->cols < 0) {
  103. G_fatal_error(_("Invalid coordinates: negative number of columns"));
  104. }
  105. if (cellhd->rows < 0) {
  106. G_fatal_error(_("Invalid coordinates: negative number of rows"));
  107. }
  108. /* (re)compute the resolutions */
  109. old_res = cellhd->ns_res;
  110. cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
  111. if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
  112. G_verbose_message(_("NS resolution has been changed"));
  113. old_res = cellhd->ew_res;
  114. cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
  115. if (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
  116. G_verbose_message(_("EW resolution has been changed"));
  117. if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
  118. G_verbose_message(_("NS and EW resolutions are different"));
  119. ll_check_ns(cellhd);
  120. ll_check_ew(cellhd);
  121. }
  122. /*!
  123. * \brief Adjust cell header for 3D values.
  124. *
  125. * This function fills in missing parts of the input cell header (or
  126. * region). It also makes projection-specific adjustments. The
  127. * <i>cellhd</i> structure must have its <i>north, south, east,
  128. * west</i>, and <i>proj</i> fields set.
  129. *
  130. * If <i>row_flag</i> is true, then the north-south resolution is computed
  131. * from the number of <i>rows</i> in the <i>cellhd</i> structure.
  132. * Otherwise the number of <i>rows</i> is computed from the north-south
  133. * resolution in the structure, similarly for <i>col_flag</i> and the
  134. * number of columns and the east-west resolution.
  135. *
  136. * If <i>depth_flag</i> is true, top-bottom resolution is calculated
  137. * from depths.
  138. * If <i>depth_flag</i> are false, number of depths is calculated from
  139. * top-bottom resolution.
  140. *
  141. * \warning This function can cause segmentation fault without any warning
  142. * when it is called with Cell_head top and bottom set to zero.
  143. *
  144. * \param[in,out] cellhd pointer to Cell_head structure
  145. * \param row_flag compute n-s resolution
  146. * \param col_flag compute e-w resolution
  147. * \param depth_flag compute t-b resolution
  148. */
  149. void G_adjust_Cell_head3(struct Cell_head *cellhd, int row_flag,
  150. int col_flag, int depth_flag)
  151. {
  152. double old_res;
  153. if (!row_flag) {
  154. if (cellhd->ns_res <= 0)
  155. G_fatal_error(_("Illegal n-s resolution value: %g"),
  156. cellhd->ns_res);
  157. if (cellhd->ns_res3 <= 0)
  158. G_fatal_error(_("Illegal n-s resolution value for 3D: %g"),
  159. cellhd->ns_res3);
  160. }
  161. else {
  162. if (cellhd->rows <= 0)
  163. G_fatal_error(_("Illegal number of rows: %d"
  164. " (resolution is %g)"),
  165. cellhd->rows, cellhd->ns_res);
  166. if (cellhd->rows3 <= 0)
  167. G_fatal_error(_("Illegal number of rows for 3D: %d"
  168. " (resolution is %g)"),
  169. cellhd->rows3, cellhd->ns_res3);
  170. }
  171. if (!col_flag) {
  172. if (cellhd->ew_res <= 0)
  173. G_fatal_error(_("Illegal e-w resolution value: %g"),
  174. cellhd->ew_res);
  175. if (cellhd->ew_res3 <= 0)
  176. G_fatal_error(_("Illegal e-w resolution value for 3D: %g"),
  177. cellhd->ew_res3);
  178. }
  179. else {
  180. if (cellhd->cols <= 0)
  181. G_fatal_error(_("Illegal number of columns: %d"
  182. " (resolution is %g)"),
  183. cellhd->cols, cellhd->ew_res);
  184. if (cellhd->cols3 <= 0)
  185. G_fatal_error(_("Illegal number of columns for 3D: %d"
  186. " (resolution is %g)"),
  187. cellhd->cols3, cellhd->ew_res3);
  188. }
  189. if (!depth_flag) {
  190. if (cellhd->tb_res <= 0)
  191. G_fatal_error(_("Illegal t-b resolution value: %g"),
  192. cellhd->tb_res);
  193. }
  194. else {
  195. if (cellhd->depths <= 0)
  196. G_fatal_error(_("Illegal depths value: %d"), cellhd->depths);
  197. }
  198. /* check the edge values */
  199. if (cellhd->north <= cellhd->south) {
  200. if (cellhd->proj == PROJECTION_LL)
  201. G_fatal_error(_("North must be north of South,"
  202. " but %g (north) <= %g (south"),
  203. cellhd->north, cellhd->south);
  204. else
  205. G_fatal_error(_("North must be larger than South,"
  206. " but %g (north) <= %g (south"),
  207. cellhd->north, cellhd->south);
  208. }
  209. ll_wrap(cellhd);
  210. if (cellhd->east <= cellhd->west)
  211. G_fatal_error(_("East must be larger than West,"
  212. " but %g (east) <= %g (west)"),
  213. cellhd->east, cellhd->west);
  214. if (cellhd->top <= cellhd->bottom)
  215. G_fatal_error(_("Top must be larger than Bottom,"
  216. " but %g (top) <= %g (bottom)"),
  217. cellhd->top, cellhd->bottom);
  218. /* compute rows and columns, if not set */
  219. if (!row_flag) {
  220. cellhd->rows =
  221. (cellhd->north - cellhd->south +
  222. cellhd->ns_res / 2.0) / cellhd->ns_res;
  223. if (cellhd->rows == 0)
  224. cellhd->rows = 1;
  225. cellhd->rows3 =
  226. (cellhd->north - cellhd->south +
  227. cellhd->ns_res3 / 2.0) / cellhd->ns_res3;
  228. if (cellhd->rows3 == 0)
  229. cellhd->rows3 = 1;
  230. }
  231. if (!col_flag) {
  232. cellhd->cols =
  233. (cellhd->east - cellhd->west +
  234. cellhd->ew_res / 2.0) / cellhd->ew_res;
  235. if (cellhd->cols == 0)
  236. cellhd->cols = 1;
  237. cellhd->cols3 =
  238. (cellhd->east - cellhd->west +
  239. cellhd->ew_res3 / 2.0) / cellhd->ew_res3;
  240. if (cellhd->cols3 == 0)
  241. cellhd->cols3 = 1;
  242. }
  243. if (!depth_flag) {
  244. cellhd->depths =
  245. (cellhd->top - cellhd->bottom +
  246. cellhd->tb_res / 2.0) / cellhd->tb_res;
  247. if (cellhd->depths == 0)
  248. cellhd->depths = 1;
  249. }
  250. if (cellhd->cols < 0 || cellhd->cols3 < 0) {
  251. G_fatal_error(_("Invalid coordinates: negative number of columns"));
  252. }
  253. if (cellhd->rows < 0 || cellhd->rows3 < 0) {
  254. G_fatal_error(_("Invalid coordinates: negative number of rows"));
  255. }
  256. if (cellhd->depths < 0) {
  257. G_fatal_error(_("Invalid coordinates: negative number of depths"));
  258. }
  259. /* (re)compute the resolutions */
  260. old_res = cellhd->ns_res;
  261. cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
  262. if (old_res > 0 && fabs(old_res - cellhd->ns_res) / old_res > 0.01)
  263. G_verbose_message(_("NS resolution has been changed"));
  264. old_res = cellhd->ew_res;
  265. cellhd->ew_res = (cellhd->east - cellhd->west) / cellhd->cols;
  266. if (old_res > 0 && fabs(old_res - cellhd->ew_res) / old_res > 0.01)
  267. G_verbose_message(_("EW resolution has been changed"));
  268. if (fabs(cellhd->ns_res - cellhd->ew_res) / cellhd->ns_res > 0.01)
  269. G_verbose_message(_("NS and EW resolutions are different"));
  270. ll_check_ns(cellhd);
  271. ll_check_ew(cellhd);
  272. cellhd->ns_res3 = (cellhd->north - cellhd->south) / cellhd->rows3;
  273. cellhd->ew_res3 = (cellhd->east - cellhd->west) / cellhd->cols3;
  274. cellhd->tb_res = (cellhd->top - cellhd->bottom) / cellhd->depths;
  275. }
  276. static int ll_wrap(struct Cell_head *cellhd)
  277. {
  278. double shift;
  279. /* for lat/lon, force east larger than west, try to wrap to -180, 180 */
  280. if (cellhd->proj != PROJECTION_LL)
  281. return 0;
  282. if (cellhd->east <= cellhd->west) {
  283. G_warning(_("East (%.15g) is not larger than West (%.15g)"),
  284. cellhd->east, cellhd->west);
  285. while (cellhd->east <= cellhd->west)
  286. cellhd->east += 360.0;
  287. }
  288. /* with east larger than west,
  289. * any 360 degree W-E extent can be represented within -360, 360
  290. * but not within -180, 180 */
  291. /* try to shift to within -180, 180 */
  292. shift = 0;
  293. while (cellhd->west + shift >= 180) {
  294. shift -= 360.0;
  295. }
  296. while (cellhd->east + shift <= -180) {
  297. shift += 360.0;
  298. }
  299. /* try to shift to within -360, 360 */
  300. while (cellhd->east + shift > 360) {
  301. shift -= 360.0;
  302. }
  303. while (cellhd->west + shift <= -360) {
  304. shift += 360.0;
  305. }
  306. if (shift) {
  307. cellhd->west += shift;
  308. cellhd->east += shift;
  309. }
  310. /* very liberal thresholds */
  311. if (cellhd->north > 90.0 + LL_TOLERANCE)
  312. G_fatal_error(_("Illegal latitude for North: %g"), cellhd->north);
  313. if (cellhd->south < -90.0 - LL_TOLERANCE)
  314. G_fatal_error(_("Illegal latitude for South: %g"), cellhd->south);
  315. #if 0
  316. /* disabled: allow W-E extents larger than 360 degree e.g. for display */
  317. if (cellhd->west < -360.0 - LL_TOLERANCE) {
  318. G_debug(1, "East: %g", cellhd->east);
  319. G_fatal_error(_("Illegal longitude for West: %g"), cellhd->west);
  320. }
  321. if (cellhd->east > 360.0 + LL_TOLERANCE) {
  322. G_debug(1, "West: %g", cellhd->west);
  323. G_fatal_error(_("Illegal longitude for East: %g"), cellhd->east);
  324. }
  325. #endif
  326. return 1;
  327. }
  328. static int ll_check_ns(struct Cell_head *cellhd)
  329. {
  330. int lladjust;
  331. double diff;
  332. int ncells;
  333. /* lat/lon checks */
  334. if (cellhd->proj != PROJECTION_LL)
  335. return 0;
  336. lladjust = 0;
  337. G_debug(3, "ll_check_ns: epsilon: %g", llepsilon);
  338. /* North, South: allow a half cell spill-over */
  339. diff = (cellhd->north - cellhd->south) / cellhd->ns_res;
  340. ncells = (int) (diff + 0.5);
  341. diff -= ncells;
  342. if ((diff < 0 && diff < -fpepsilon) ||
  343. (diff > 0 && diff > fpepsilon)) {
  344. G_verbose_message(_("NS extent does not match NS resolution: %g cells difference"),
  345. diff);
  346. }
  347. /* north */
  348. diff = (cellhd->north - 90) / cellhd->ns_res;
  349. if (diff < 0)
  350. diff = -diff;
  351. if (cellhd->north < 90.0 && diff < 1.0 ) {
  352. G_verbose_message(_("%g cells missing to reach 90 degree north"),
  353. diff);
  354. if (diff < llepsilon && diff > fpepsilon) {
  355. G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
  356. cellhd->north - 90.0);
  357. /* check only, do not modify
  358. cellhd->north = 90.0;
  359. lladjust = 1;
  360. */
  361. }
  362. }
  363. if (cellhd->north > 90.0) {
  364. if (diff <= 0.5 + llepsilon) {
  365. G_important_message(_("90 degree north is exceeded by %g cells"),
  366. diff);
  367. if (diff < llepsilon && diff > fpepsilon) {
  368. G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
  369. cellhd->north - 90.0);
  370. G_debug(1, "North of north in seconds: %g",
  371. (cellhd->north - 90.0) * 3600);
  372. /* check only, do not modify
  373. cellhd->north = 90.0;
  374. lladjust = 1;
  375. */
  376. }
  377. diff = diff - 0.5;
  378. if (diff < 0)
  379. diff = -diff;
  380. if (diff < llepsilon && diff > fpepsilon) {
  381. G_verbose_message(_("Subtle input data rounding error of north boundary (%g)"),
  382. cellhd->north - 90.0 - cellhd->ns_res / 2.0);
  383. G_debug(1, "North of north + 0.5 cells in seconds: %g",
  384. (cellhd->north - 90.0 - cellhd->ns_res / 2.0) * 3600);
  385. /* check only, do not modify
  386. cellhd->north = 90.0 + cellhd->ns_res / 2.0;
  387. lladjust = 1;
  388. */
  389. }
  390. }
  391. else
  392. G_fatal_error(_("Illegal latitude for North: %g"), cellhd->north);
  393. }
  394. /* south */
  395. diff = (cellhd->south + 90) / cellhd->ns_res;
  396. if (diff < 0)
  397. diff = -diff;
  398. if (cellhd->south > -90.0 && diff < 1.0 ) {
  399. G_verbose_message(_("%g cells missing to reach 90 degree south"),
  400. diff);
  401. if (diff < llepsilon && diff > fpepsilon) {
  402. G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
  403. cellhd->south + 90.0);
  404. /* check only, do not modify
  405. cellhd->south = -90.0;
  406. lladjust = 1;
  407. */
  408. }
  409. }
  410. if (cellhd->south < -90.0) {
  411. if (diff <= 0.5 + llepsilon) {
  412. G_important_message(_("90 degree south is exceeded by %g cells"),
  413. diff);
  414. if (diff < llepsilon && diff > fpepsilon) {
  415. G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
  416. cellhd->south + 90);
  417. G_debug(1, "South of south in seconds: %g",
  418. (-cellhd->south - 90) * 3600);
  419. /* check only, do not modify
  420. cellhd->south = -90.0;
  421. lladjust = 1;
  422. */
  423. }
  424. diff = diff - 0.5;
  425. if (diff < 0)
  426. diff = -diff;
  427. if (diff < llepsilon && diff > fpepsilon) {
  428. G_verbose_message(_("Subtle input data rounding error of south boundary (%g)"),
  429. cellhd->south + 90 + cellhd->ns_res / 2.0);
  430. G_debug(1, "South of south + 0.5 cells in seconds: %g",
  431. (-cellhd->south - 90 - cellhd->ns_res / 2.0) * 3600);
  432. /* check only, do not modify
  433. cellhd->south = -90.0 - cellhd->ns_res / 2.0;
  434. lladjust = 1;
  435. */
  436. }
  437. }
  438. else
  439. G_fatal_error(_("Illegal latitude for South: %g"), cellhd->south);
  440. }
  441. if (lladjust)
  442. cellhd->ns_res = (cellhd->north - cellhd->south) / cellhd->rows;
  443. return lladjust;
  444. }
  445. static int ll_check_ew(struct Cell_head *cellhd)
  446. {
  447. int lladjust;
  448. double diff;
  449. int ncells;
  450. /* lat/lon checks */
  451. if (cellhd->proj != PROJECTION_LL)
  452. return 0;
  453. lladjust = 0;
  454. G_debug(3, "ll_check_ew: epsilon: %g", llepsilon);
  455. /* west - east, no adjustment */
  456. diff = (cellhd->east - cellhd->west) / cellhd->ew_res;
  457. ncells = (int) (diff + 0.5);
  458. diff -= ncells;
  459. if ((diff < 0 && diff < -fpepsilon) ||
  460. (diff > 0 && diff > fpepsilon)) {
  461. G_verbose_message(_("EW extent does not match EW resolution: %g cells difference"),
  462. diff);
  463. }
  464. if (cellhd->east - cellhd->west > 360.0) {
  465. diff = (cellhd->east - cellhd->west - 360.0) / cellhd->ew_res;
  466. if (diff > fpepsilon)
  467. G_important_message(_("360 degree EW extent is exceeded by %g cells"
  468. " (East: %g, West: %g)"),
  469. diff, cellhd->east, cellhd->west);
  470. }
  471. else if (cellhd->east - cellhd->west < 360.0) {
  472. diff = (360.0 - (cellhd->east - cellhd->west)) / cellhd->ew_res;
  473. if (diff < 1.0 && diff > fpepsilon)
  474. G_verbose_message(_("%g cells missing to cover 360 degree EW extent"),
  475. diff);
  476. }
  477. return lladjust;
  478. }
  479. /*!
  480. * \brief Adjust window for lat/lon.
  481. *
  482. * This function tries to automatically fix fp precision issues and
  483. * adjust rounding errors for lat/lon.
  484. *
  485. * <b>Note:</b> 3D values are not adjusted.
  486. *
  487. * \param[in,out] cellhd pointer to Cell_head structure
  488. * \return 1 if window was adjusted
  489. * \return 0 if window was not adjusted
  490. */
  491. int G_adjust_window_ll(struct Cell_head *cellhd)
  492. {
  493. int ll_adjust, res_adj;
  494. double dsec, dsec2;
  495. char buf[100], buf2[100];
  496. double diff;
  497. double old, new;
  498. struct Cell_head cellhds; /* everything in seconds, not degrees */
  499. /* lat/lon checks */
  500. if (cellhd->proj != PROJECTION_LL)
  501. return 0;
  502. /* put everything through ll_format + ll_scan */
  503. G_llres_format(cellhd->ns_res, buf);
  504. if (G_llres_scan(buf, &new) != 1)
  505. G_fatal_error(_("Invalid NS resolution"));
  506. cellhd->ns_res = new;
  507. G_llres_format(cellhd->ew_res, buf);
  508. if (G_llres_scan(buf, &new) != 1)
  509. G_fatal_error(_("Invalid EW resolution"));
  510. cellhd->ew_res = new;
  511. G_lat_format(cellhd->north, buf);
  512. if (G_lat_scan(buf, &new) != 1)
  513. G_fatal_error(_("Invalid North"));
  514. cellhd->north = new;
  515. G_lat_format(cellhd->south, buf);
  516. if (G_lat_scan(buf, &new) != 1)
  517. G_fatal_error(_("Invalid South"));
  518. cellhd->south = new;
  519. G_lon_format(cellhd->west, buf);
  520. if (G_lon_scan(buf, &new) != 1)
  521. G_fatal_error(_("Invalid West"));
  522. cellhd->west = new;
  523. G_lon_format(cellhd->east, buf);
  524. if (G_lon_scan(buf, &new) != 1)
  525. G_fatal_error(_("Invalid East"));
  526. cellhd->east = new;
  527. /* convert to seconds */
  528. cellhds = *cellhd;
  529. old = cellhds.ns_res * 3600;
  530. sprintf(buf, "%f", old);
  531. sscanf(buf, "%lf", &new);
  532. cellhds.ns_res = new;
  533. old = cellhds.ew_res * 3600;
  534. sprintf(buf, "%f", old);
  535. sscanf(buf, "%lf", &new);
  536. cellhds.ew_res = new;
  537. old = cellhds.north * 3600;
  538. sprintf(buf, "%f", old);
  539. sscanf(buf, "%lf", &new);
  540. cellhds.north = new;
  541. old = cellhds.south * 3600;
  542. sprintf(buf, "%f", old);
  543. sscanf(buf, "%lf", &new);
  544. cellhds.south = new;
  545. old = cellhds.west * 3600;
  546. sprintf(buf, "%f", old);
  547. sscanf(buf, "%lf", &new);
  548. cellhds.west = new;
  549. old = cellhds.east * 3600;
  550. sprintf(buf, "%f", old);
  551. sscanf(buf, "%lf", &new);
  552. cellhds.east = new;
  553. ll_adjust = 0;
  554. /* N - S */
  555. /* resolution */
  556. res_adj = 0;
  557. old = cellhds.ns_res;
  558. if (old > 0.4) {
  559. /* round to nearest 0.1 sec */
  560. dsec = old * 10;
  561. dsec2 = floor(dsec + 0.5);
  562. new = dsec2 / 10;
  563. diff = fabs(dsec2 - dsec) / dsec;
  564. if (diff > 0 && diff < llepsilon) {
  565. G_llres_format(old / 3600, buf);
  566. G_llres_format(new / 3600, buf2);
  567. if (strcmp(buf, buf2))
  568. G_verbose_message(_("NS resolution rounded from %s to %s"),
  569. buf, buf2);
  570. ll_adjust = 1;
  571. res_adj = 1;
  572. cellhds.ns_res = new;
  573. }
  574. }
  575. if (res_adj) {
  576. double n_off, s_off;
  577. old = cellhds.north;
  578. dsec = old * 10;
  579. dsec2 = floor(dsec + 0.5);
  580. diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
  581. n_off = diff;
  582. old = cellhds.south;
  583. dsec = old * 10;
  584. dsec2 = floor(dsec + 0.5);
  585. diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
  586. s_off = diff;
  587. if (n_off < llepsilon || n_off <= s_off) {
  588. old = cellhds.north;
  589. dsec = old * 10;
  590. dsec2 = floor(dsec + 0.5);
  591. new = dsec2 / 10;
  592. diff = n_off;
  593. if (diff > 0 && diff < llepsilon) {
  594. G_lat_format(old / 3600, buf);
  595. G_lat_format(new / 3600, buf2);
  596. if (strcmp(buf, buf2))
  597. G_verbose_message(_("North rounded from %s to %s"),
  598. buf, buf2);
  599. cellhds.north = new;
  600. }
  601. old = cellhds.south;
  602. new = cellhds.north - cellhds.ns_res * cellhds.rows;
  603. diff = fabs(new - old) / cellhds.ns_res;
  604. if (diff > 0) {
  605. G_lat_format(old / 3600, buf);
  606. G_lat_format(new / 3600, buf2);
  607. if (strcmp(buf, buf2))
  608. G_verbose_message(_("South adjusted from %s to %s"),
  609. buf, buf2);
  610. }
  611. cellhds.south = new;
  612. }
  613. else {
  614. old = cellhds.south;
  615. dsec = old * 10;
  616. dsec2 = floor(dsec + 0.5);
  617. new = dsec2 / 10;
  618. diff = s_off;
  619. if (diff > 0 && diff < llepsilon) {
  620. G_lat_format(old / 3600, buf);
  621. G_lat_format(new / 3600, buf2);
  622. if (strcmp(buf, buf2))
  623. G_verbose_message(_("South rounded from %s to %s"),
  624. buf, buf2);
  625. cellhds.south = new;
  626. }
  627. old = cellhds.north;
  628. new = cellhds.south + cellhds.ns_res * cellhds.rows;
  629. diff = fabs(new - old) / cellhds.ns_res;
  630. if (diff > 0) {
  631. G_lat_format(old / 3600, buf);
  632. G_lat_format(new / 3600, buf2);
  633. if (strcmp(buf, buf2))
  634. G_verbose_message(_("North adjusted from %s to %s"),
  635. buf, buf2);
  636. }
  637. cellhds.north = new;
  638. }
  639. }
  640. else {
  641. old = cellhds.north;
  642. dsec = old * 10;
  643. dsec2 = floor(dsec + 0.5);
  644. new = dsec2 / 10;
  645. diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
  646. if (diff > 0 && diff < llepsilon) {
  647. G_lat_format(old / 3600, buf);
  648. G_lat_format(new / 3600, buf2);
  649. if (strcmp(buf, buf2))
  650. G_verbose_message(_("North rounded from %s to %s"),
  651. buf, buf2);
  652. ll_adjust = 1;
  653. cellhds.north = new;
  654. }
  655. old = cellhds.south;
  656. dsec = old * 10;
  657. dsec2 = floor(dsec + 0.5);
  658. new = dsec2 / 10;
  659. diff = fabs(dsec2 - dsec) / (cellhds.ns_res * 10);
  660. if (diff > 0 && diff < llepsilon) {
  661. G_lat_format(old / 3600, buf);
  662. G_lat_format(new / 3600, buf2);
  663. if (strcmp(buf, buf2))
  664. G_verbose_message(_("South rounded from %s to %s"),
  665. buf, buf2);
  666. ll_adjust = 1;
  667. cellhds.south = new;
  668. }
  669. }
  670. cellhds.ns_res = (cellhds.north - cellhds.south) / cellhds.rows;
  671. /* E - W */
  672. /* resolution */
  673. res_adj = 0;
  674. old = cellhds.ew_res;
  675. if (old > 0.4) {
  676. /* round to nearest 0.1 sec */
  677. dsec = old * 10;
  678. dsec2 = floor(dsec + 0.5);
  679. new = dsec2 / 10;
  680. diff = fabs(dsec2 - dsec) / dsec;
  681. if (diff > 0 && diff < llepsilon) {
  682. G_llres_format(old / 3600, buf);
  683. G_llres_format(new / 3600, buf2);
  684. if (strcmp(buf, buf2))
  685. G_verbose_message(_("EW resolution rounded from %s to %s"),
  686. buf, buf2);
  687. ll_adjust = 1;
  688. res_adj = 1;
  689. cellhds.ew_res = new;
  690. }
  691. }
  692. if (res_adj) {
  693. double w_off, e_off;
  694. old = cellhds.west;
  695. dsec = old * 10;
  696. dsec2 = floor(dsec + 0.5);
  697. diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
  698. w_off = diff;
  699. old = cellhds.east;
  700. dsec = old * 10;
  701. dsec2 = floor(dsec + 0.5);
  702. diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
  703. e_off = diff;
  704. if (w_off < llepsilon || w_off <= e_off) {
  705. old = cellhds.west;
  706. dsec = old * 10;
  707. dsec2 = floor(dsec + 0.5);
  708. new = dsec2 / 10;
  709. diff = w_off;
  710. if (diff > 0 && diff < llepsilon) {
  711. G_lon_format(old / 3600, buf);
  712. G_lon_format(new / 3600, buf2);
  713. if (strcmp(buf, buf2))
  714. G_verbose_message(_("West rounded from %s to %s"),
  715. buf, buf2);
  716. cellhds.west = new;
  717. }
  718. old = cellhds.east;
  719. new = cellhds.west + cellhds.ew_res * cellhds.cols;
  720. diff = fabs(new - old) / cellhds.ew_res;
  721. if (diff > 0) {
  722. G_lon_format(old / 3600, buf);
  723. G_lon_format(new / 3600, buf2);
  724. if (strcmp(buf, buf2))
  725. G_verbose_message(_("East adjusted from %s to %s"),
  726. buf, buf2);
  727. }
  728. cellhds.east = new;
  729. }
  730. else {
  731. old = cellhds.east;
  732. dsec = old * 10;
  733. dsec2 = floor(dsec + 0.5);
  734. new = dsec2 / 10;
  735. diff = e_off;
  736. if (diff > 0 && diff < llepsilon) {
  737. G_lon_format(old / 3600, buf);
  738. G_lon_format(new / 3600, buf2);
  739. if (strcmp(buf, buf2))
  740. G_verbose_message(_("East rounded from %s to %s"),
  741. buf, buf2);
  742. cellhds.east = new;
  743. }
  744. old = cellhds.west;
  745. new = cellhds.east - cellhds.ew_res * cellhds.cols;
  746. diff = fabs(new - cellhds.west) / cellhds.ew_res;
  747. if (diff > 0) {
  748. G_lon_format(old / 3600, buf);
  749. G_lon_format(new / 3600, buf2);
  750. if (strcmp(buf, buf2))
  751. G_verbose_message(_("West adjusted from %s to %s"),
  752. buf, buf2);
  753. }
  754. cellhds.west = new;
  755. }
  756. }
  757. else {
  758. old = cellhds.west;
  759. dsec = old * 10;
  760. dsec2 = floor(dsec + 0.5);
  761. new = dsec2 / 10;
  762. diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
  763. if (diff > 0 && diff < llepsilon) {
  764. G_lon_format(old / 3600, buf);
  765. G_lon_format(new / 3600, buf2);
  766. if (strcmp(buf, buf2))
  767. G_verbose_message(_("West rounded from %s to %s"),
  768. buf, buf2);
  769. ll_adjust = 1;
  770. cellhds.west = new;
  771. }
  772. old = cellhds.east;
  773. dsec = old * 10;
  774. dsec2 = floor(dsec + 0.5);
  775. new = dsec2 / 10;
  776. diff = fabs(dsec2 - dsec) / (cellhds.ew_res * 10);
  777. if (diff > 0 && diff < llepsilon) {
  778. G_lon_format(old / 3600, buf);
  779. G_lon_format(new / 3600, buf2);
  780. if (strcmp(buf, buf2))
  781. G_verbose_message(_("East rounded from %s to %s"),
  782. buf, buf2);
  783. ll_adjust = 1;
  784. cellhds.east = new;
  785. }
  786. }
  787. cellhds.ew_res = (cellhds.east - cellhds.west) / cellhds.cols;
  788. cellhd->ns_res = cellhds.ns_res / 3600;
  789. cellhd->ew_res = cellhds.ew_res / 3600;
  790. cellhd->north = cellhds.north / 3600;
  791. cellhd->south = cellhds.south / 3600;
  792. cellhd->west = cellhds.west / 3600;
  793. cellhd->east = cellhds.east / 3600;
  794. return ll_adjust;
  795. }