cnversions.c 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472
  1. #include <math.h>
  2. #include <grass/gis.h>
  3. #include <grass/display.h>
  4. /****** OLD CODE
  5. * #include "windround.h"
  6. **********/
  7. /* D_do_conversions(window, t, b, l, r)
  8. * struct Cell_head *window ;
  9. * int t, b, l, r ;
  10. *
  11. * Sets up conversion coefficients to translate between three
  12. * coordinate systems:
  13. *
  14. * 1. Screen coordinates (given by t, b, l, r values)
  15. * 2. UTM coordinates (given by values in window structure)
  16. * 3. Window array coors (given by values in window structure)
  17. *
  18. * Once D_do_conversions is called, lots of conversion coefficients
  19. * and conversion routines are available.
  20. *
  21. * Calls to convert row and column (x and y) values in one system to
  22. * another system are available. In addition calls which return the
  23. * conversion coefficients are alos provided.
  24. */
  25. struct vector
  26. {
  27. double x, y;
  28. };
  29. struct rectangle
  30. {
  31. double west;
  32. double east;
  33. double south;
  34. double north;
  35. struct vector size;
  36. };
  37. /* Bounding rectangles */
  38. static struct rectangle D; /* Display coordinates, pixels, (0,0) towards NW */
  39. static struct rectangle A; /* Map array coordinates, integers, (0,0) towards NW */
  40. static struct rectangle U; /* UTM coordinates, meters, (0,0) towards SW */
  41. /* Conversion factors */
  42. static struct vector D_to_A_conv; /* Display to Array */
  43. static struct vector A_to_U_conv; /* Array to UTM */
  44. static struct vector U_to_D_conv; /* UTM to Display */
  45. /* others */
  46. static int is_lat_lon;
  47. /*!
  48. * \brief initialize conversions
  49. *
  50. * The relationship between the
  51. * earth <b>region</b> and the <b>top, bottom, left</b>, and <b>right</b>
  52. * screen coordinates is established, which then allows conversions between all
  53. * three coordinate systems to be performed.
  54. * Note this routine is called by <i>D_setup.</i>
  55. *
  56. * \param region
  57. * \param top
  58. * \param bottom
  59. * \param left
  60. * \param right
  61. * \return int
  62. */
  63. static void calc_size(struct rectangle *rect)
  64. {
  65. rect->size.x = rect->east - rect->west;
  66. rect->size.y = rect->south - rect->north;
  67. }
  68. static void calc_conv(struct vector *conv,
  69. const struct vector *src, const struct vector *dst)
  70. {
  71. conv->x = dst->x / src->x;
  72. conv->y = dst->y / src->y;
  73. }
  74. static void fit_aspect(struct rectangle *rect, const struct rectangle *ref)
  75. {
  76. struct vector conv;
  77. double scale, size, delta;
  78. calc_conv(&conv, &rect->size, &ref->size);
  79. if (fabs(conv.y) > fabs(conv.x)) {
  80. scale = fabs(conv.y) / fabs(conv.x);
  81. size = rect->size.x / scale;
  82. delta = rect->size.x - size;
  83. rect->west += delta/2;
  84. rect->east -= delta/2;
  85. rect->size.x = size;
  86. }
  87. else {
  88. scale = fabs(conv.x) / fabs(conv.y);
  89. size = rect->size.y / scale;
  90. delta = rect->size.y - size;
  91. rect->north += delta/2;
  92. rect->south -= delta/2;
  93. rect->size.y = size;
  94. }
  95. }
  96. void D_update_conversions(void)
  97. {
  98. calc_conv(&D_to_A_conv, &D.size, &A.size);
  99. calc_conv(&A_to_U_conv, &A.size, &U.size);
  100. calc_conv(&U_to_D_conv, &U.size, &D.size);
  101. }
  102. void D_fit_d_to_u(void)
  103. {
  104. fit_aspect(&D, &U);
  105. }
  106. void D_fit_u_to_d(void)
  107. {
  108. fit_aspect(&U, &D);
  109. }
  110. void D_show_conversions(void)
  111. {
  112. fprintf(stderr,
  113. " D_w %10.1f D_e %10.1f D_s %10.1f D_n %10.1f\n",
  114. D.west, D.east, D.south, D.north);
  115. fprintf(stderr,
  116. " A_w %10.1f A_e %10.1f A_s %10.1f A_n %10.1f\n",
  117. A.west, A.east, A.south, A.north);
  118. fprintf(stderr,
  119. " U_w %10.1f U_e %10.1f U_s %10.1f U_n %10.1f\n",
  120. U.west, U.east, U.south, U.north);
  121. fprintf(stderr,
  122. " D_x %10.1f D_y %10.1f\n" "\n", D.size.x, D.size.y);
  123. fprintf(stderr,
  124. " A_x %10.1f A_y %10.1f\n" "\n", A.size.x, A.size.y);
  125. fprintf(stderr,
  126. " U_x %10.1f U_y %10.1f\n" "\n", U.size.x, U.size.y);
  127. fprintf(stderr, " D_to_A_conv.x %10.1f D_to_A_conv.y %10.1f \n",
  128. D_to_A_conv.x, D_to_A_conv.y);
  129. fprintf(stderr, " A_to_U_conv.x %10.1f A_to_U_conv.y %10.1f \n",
  130. A_to_U_conv.x, A_to_U_conv.y);
  131. fprintf(stderr, " U_to_D_conv.x %10.1f U_to_D_conv.y %10.1f \n",
  132. U_to_D_conv.x, U_to_D_conv.y);
  133. }
  134. void D_do_conversions(const struct Cell_head *window,
  135. double t, double b, double l, double r)
  136. {
  137. D_set_region(window);
  138. D_set_dst(t, b, l, r);
  139. D_fit_d_to_u();
  140. D_update_conversions();
  141. #ifdef DEBUG
  142. D_show_conversions();
  143. #endif /* DEBUG */
  144. }
  145. int D_is_lat_lon(void) { return (is_lat_lon); }
  146. double D_get_d_to_a_xconv(void) { return (D_to_A_conv.x); }
  147. double D_get_d_to_a_yconv(void) { return (D_to_A_conv.y); }
  148. double D_get_d_to_u_xconv(void) { return (1/U_to_D_conv.x); }
  149. double D_get_d_to_u_yconv(void) { return (1/U_to_D_conv.y); }
  150. double D_get_a_to_u_xconv(void) { return (A_to_U_conv.x); }
  151. double D_get_a_to_u_yconv(void) { return (A_to_U_conv.y); }
  152. double D_get_a_to_d_xconv(void) { return (1/D_to_A_conv.x); }
  153. double D_get_a_to_d_yconv(void) { return (1/D_to_A_conv.y); }
  154. double D_get_u_to_d_xconv(void) { return (U_to_D_conv.x); }
  155. double D_get_u_to_d_yconv(void) { return (U_to_D_conv.y); }
  156. double D_get_u_to_a_xconv(void) { return (1/A_to_U_conv.x); }
  157. double D_get_u_to_a_yconv(void) { return (1/A_to_U_conv.y); }
  158. double D_get_ns_resolution(void) { return D_get_a_to_u_yconv(); }
  159. double D_get_ew_resolution(void) { return D_get_a_to_u_xconv(); }
  160. double D_get_u_west(void) { return (U.west); }
  161. double D_get_u_east(void) { return (U.east); }
  162. double D_get_u_north(void) { return (U.north); }
  163. double D_get_u_south(void) { return (U.south); }
  164. double D_get_a_west(void) { return (A.west); }
  165. double D_get_a_east(void) { return (A.east); }
  166. double D_get_a_north(void) { return (A.north); }
  167. double D_get_a_south(void) { return (A.south); }
  168. double D_get_d_west(void) { return (D.west); }
  169. double D_get_d_east(void) { return (D.east); }
  170. double D_get_d_north(void) { return (D.north); }
  171. double D_get_d_south(void) { return (D.south); }
  172. void D_set_region(const struct Cell_head *window)
  173. {
  174. D_set_src(window->north, window->south, window->west, window->east);
  175. D_set_grid(0, window->rows, 0, window->cols);
  176. is_lat_lon = (window->proj == PROJECTION_LL);
  177. }
  178. void D_set_src(double t, double b, double l, double r)
  179. {
  180. U.north = t;
  181. U.south = b;
  182. U.west = l;
  183. U.east = r;
  184. calc_size(&U);
  185. }
  186. void D_get_src(double *t, double *b, double *l, double *r)
  187. {
  188. *t = U.north;
  189. *b = U.south;
  190. *l = U.west;
  191. *r = U.east;
  192. }
  193. void D_set_grid(int t, int b, int l, int r)
  194. {
  195. A.north = t;
  196. A.south = b;
  197. A.west = l;
  198. A.east = r;
  199. calc_size(&A);
  200. }
  201. void D_get_grid(int *t, int *b, int *l, int *r)
  202. {
  203. *t = A.north;
  204. *b = A.south;
  205. *l = A.west;
  206. *r = A.east;
  207. }
  208. void D_set_dst(double t, double b, double l, double r)
  209. {
  210. D.north = t;
  211. D.south = b;
  212. D.west = l;
  213. D.east = r;
  214. calc_size(&D);
  215. }
  216. void D_get_dst(double *t, double *b, double *l, double *r)
  217. {
  218. *t = D.north;
  219. *b = D.south;
  220. *l = D.west;
  221. *r = D.east;
  222. }
  223. void D_get_u(double x[2][2])
  224. {
  225. x[0][0] = U.west;
  226. x[0][1] = U.east;
  227. x[1][0] = U.north;
  228. x[1][1] = U.south;
  229. }
  230. void D_get_a(int x[2][2])
  231. {
  232. x[0][0] = (int)A.west;
  233. x[0][1] = (int)A.east;
  234. x[1][0] = (int)A.north;
  235. x[1][1] = (int)A.south;
  236. }
  237. void D_get_d(double x[2][2])
  238. {
  239. x[0][0] = D.west;
  240. x[0][1] = D.east;
  241. x[1][0] = D.north;
  242. x[1][1] = D.south;
  243. }
  244. /*!
  245. * \brief screen to array (y)
  246. *
  247. * Returns a <i>row</i> value in the array coordinate system when provided the
  248. * corresponding <b>y</b> value in the screen coordinate system.
  249. *
  250. * \param y
  251. * \return double
  252. */
  253. double D_d_to_a_row(double D_row)
  254. {
  255. return A.north + (D_row - D.north) * D_to_A_conv.y;
  256. }
  257. /*!
  258. * \brief screen to array (x)
  259. *
  260. * Returns a <i>column</i> value in the array coordinate system when provided the
  261. * corresponding <b>x</b> value in the screen coordinate system.
  262. *
  263. * \param x
  264. * \return double
  265. */
  266. double D_d_to_a_col(double D_col)
  267. {
  268. return A.west + (D_col - D.west) * D_to_A_conv.x;
  269. }
  270. /*!
  271. * \brief screen to earth (y)
  272. *
  273. * Returns a <i>north</i> value in the earth coordinate system when provided the
  274. * corresponding <b>y</b> value in the screen coordinate system.
  275. *
  276. * \param y
  277. * \return double
  278. */
  279. double D_d_to_u_row(double D_row)
  280. {
  281. return U.north + (D_row - D.north) / U_to_D_conv.y;
  282. }
  283. /*!
  284. * \brief screen to earth (x)
  285. *
  286. * Returns an <i>east</i> value in the earth coordinate system when provided the
  287. * corresponding <b>x</b> value in the screen coordinate system.
  288. *
  289. * \param x
  290. * \return double
  291. */
  292. double D_d_to_u_col(double D_col)
  293. {
  294. return U.west + (D_col - D.west) / U_to_D_conv.x;
  295. }
  296. /*!
  297. * \brief array to earth (row)
  298. *
  299. * Returns a <i>y</i> value in the earth coordinate system when provided the
  300. * corresponding <b>row</b> value in the array coordinate system.
  301. *
  302. * \param row
  303. * \return double
  304. */
  305. double D_a_to_u_row(double A_row)
  306. {
  307. return U.north + (A_row - A.north) * A_to_U_conv.y;
  308. }
  309. /*!
  310. * \brief array to earth (column)
  311. *
  312. * Returns an <i>x</i> value in the earth coordinate system when
  313. * provided the corresponding <b>column</b> value in the array coordinate
  314. * system.
  315. *
  316. * \param column
  317. * \return double
  318. */
  319. double D_a_to_u_col(double A_col)
  320. {
  321. return U.west + (A_col - A.west) * A_to_U_conv.x;
  322. }
  323. /*!
  324. * \brief array to screen (row)
  325. *
  326. * Returns a <i>y</i> value in the screen coordinate system when provided the
  327. * corresponding <b>row</b> value in the array coordinate system.
  328. *
  329. * \param row
  330. * \return double
  331. */
  332. double D_a_to_d_row(double A_row)
  333. {
  334. return D.north + (A_row - A.north) / D_to_A_conv.y;
  335. }
  336. /*!
  337. * \brief array to screen (column)
  338. *
  339. * Returns an <i>x</i> value in the screen coordinate system when
  340. * provided the corresponding <b>column</b> value in the array coordinate
  341. * system.
  342. *
  343. * \param column
  344. * \return double
  345. */
  346. double D_a_to_d_col(double A_col)
  347. {
  348. return D.west + (A_col - A.west) / D_to_A_conv.x;
  349. }
  350. /*!
  351. * \brief earth to screen (north)
  352. *
  353. * Returns a <i>y</i> value in the screen coordinate system when provided the
  354. * corresponding <b>north</b> value in the earth coordinate system.
  355. *
  356. * \param north
  357. * \return double
  358. */
  359. double D_u_to_d_row(double U_row)
  360. {
  361. return D.north + (U_row - U.north) * U_to_D_conv.y;
  362. }
  363. /*!
  364. * \brief earth to screen (east)
  365. *
  366. * Returns an <i>x</i> value in the screen coordinate system when provided the
  367. * corresponding <b>east</b> value in the earth coordinate system.
  368. *
  369. * \param east
  370. * \return double
  371. */
  372. double D_u_to_d_col(double U_col)
  373. {
  374. return D.west + (U_col - U.west) * U_to_D_conv.x;
  375. }
  376. /*!
  377. * \brief earth to array (north)
  378. *
  379. * Returns a <i>row</i> value in the array coordinate system when provided the
  380. * corresponding <b>north</b> value in the earth coordinate system.
  381. *
  382. * \param north
  383. * \return double
  384. */
  385. double D_u_to_a_row(double U_row)
  386. {
  387. return A.north + (U_row - U.north) / A_to_U_conv.y;
  388. }
  389. /*!
  390. * \brief earth to array (east
  391. *
  392. * Returns a <i>column</i> value in the array coordinate system when provided the
  393. * corresponding <b>east</b> value in the earth coordinate system.
  394. *
  395. * \param east
  396. * \return double
  397. */
  398. double D_u_to_a_col(double U_col)
  399. {
  400. return A.west + (U_col - U.west) / A_to_U_conv.x;
  401. }