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