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