source: grass/trunk/general/g.region/printwindow.c

Last change on this file was 73852, checked in by neteler, 6 years ago

bulk fixing of typos in comments (using tools/fix_typos.sh)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id
  • Property svn:mime-type set to text/x-csrc
File size: 23.1 KB
Line 
1#include <string.h>
2#include <stdlib.h>
3#include <grass/gis.h>
4#include <grass/gprojects.h>
5#include <grass/glocale.h>
6#include "local_proto.h"
7
8#define DEG2RAD(a) ((a) * M_PI / 180.0)
9#define RAD2DEG(a) ((a) * 180.0 / M_PI)
10
11static double get_shift(double east)
12{
13 double shift;
14
15 shift = 0;
16 while (east + shift > 180)
17 shift -= 360;
18 while (east + shift < -180)
19 shift += 360;
20
21 return shift;
22}
23
24void print_window(struct Cell_head *window, int print_flag, int flat_flag)
25{
26 const char *prj, *datum, *ellps;
27 int x, width = 11;
28
29 char north[30], south[30], east[30], west[30], nsres[30], ewres[30],
30 nsres3[30], ewres3[30], tbres[30];
31 char buf[50];
32 char *sep = "\n";
33
34 double ew_dist1, ew_dist2, ns_dist1, ns_dist2;
35 double longitude, latitude;
36
37 if (print_flag & PRINT_SH)
38 x = G_projection() == PROJECTION_LL ? -1 : 0;
39 else
40 x = window->proj;
41
42 G_format_northing(window->north, north, x);
43 G_format_northing(window->south, south, x);
44 G_format_easting(window->east, east, x);
45 G_format_easting(window->west, west, x);
46 G_format_resolution(window->ew_res, ewres, x);
47 G_format_resolution(window->ew_res3, ewres3, x);
48 G_format_resolution(window->ns_res, nsres, x);
49 G_format_resolution(window->ns_res3, nsres3, x);
50 G_format_resolution(window->tb_res, tbres, -1);
51 G_begin_distance_calculations();
52
53 /* EW Dist at North edge */
54 ew_dist1 =
55 G_distance(window->east, window->north, window->west, window->north);
56 /* EW Dist at South Edge */
57 ew_dist2 =
58 G_distance(window->east, window->south, window->west, window->south);
59 /* NS Dist at East edge */
60 ns_dist1 =
61 G_distance(window->east, window->north, window->east, window->south);
62 /* NS Dist at West edge */
63 ns_dist2 =
64 G_distance(window->west, window->north, window->west, window->south);
65
66 /* width */
67 if (print_flag & PRINT_REG)
68 width = 11;
69
70 if (print_flag & PRINT_CENTER || print_flag & PRINT_MBBOX)
71 width = 16;
72
73 if (print_flag & PRINT_LL || print_flag & PRINT_NANGLE)
74 width = 18;
75
76 if (print_flag & PRINT_EXTENT)
77 width = 19;
78
79 /* flag.dist_res */
80 if (print_flag & PRINT_METERS) {
81 sprintf(ewres, "%.8f", ((ew_dist1 + ew_dist2) / 2) / window->cols);
82 G_trim_decimal(ewres);
83 sprintf(ewres3, "%.8f", ((ew_dist1 + ew_dist2) / 2) / window->cols3);
84 G_trim_decimal(ewres3);
85 sprintf(nsres, "%.8f", ((ns_dist1 + ns_dist2) / 2) / window->rows);
86 G_trim_decimal(nsres);
87 sprintf(nsres3, "%.8f", ((ns_dist1 + ns_dist2) / 2) / window->rows3);
88 G_trim_decimal(nsres3);
89 sprintf(tbres, "%.8f",
90 (window->top - window->bottom) / window->depths);
91 G_trim_decimal(tbres);
92 }
93
94 /* flag.print & flag.gprint */
95 if (print_flag & PRINT_REG) {
96 prj = G_database_projection_name();
97 if (!prj)
98 prj = "** unknown **";
99
100 if (print_flag & PRINT_SH) {
101 fprintf(stdout, "projection=%d\n", window->proj);
102 fprintf(stdout, "zone=%d\n", window->zone);
103 }
104 else {
105 fprintf(stdout, "%-*s %d (%s)\n", width, "projection:",
106 window->proj, prj);
107 fprintf(stdout, "%-*s %d\n", width, "zone:", window->zone);
108 }
109
110 /* don't print datum/ellipsoid in XY-Locations */
111 if (window->proj != 0) {
112 datum = G_database_datum_name();
113 if (!datum)
114 datum = "** unknown (default: WGS84) **";
115 ellps = G_database_ellipse_name();
116 if (!ellps)
117 ellps = "** unknown (default: WGS84) **";
118
119 /*
120 please remove before GRASS 7 released
121 backward compatibility issue
122
123 if (print_flag & PRINT_SH)
124 {
125 if (datum[0] != '*')
126 fprintf(stdout, "datum=%s\n", datum);
127 else
128 fprintf(stdout, "datum=wgs84\n");
129 if (ellps[0] != '*')
130 fprintf(stdout, "ellipsoid=%s\n", ellps);
131 else
132 fprintf(stdout, "ellipsoid=wgs84\n");
133 }
134 else
135 {
136 fprintf(stdout, "%-*s %s\n", width, "datum:", datum);
137 fprintf(stdout, "%-*s %s\n", width, "ellipsoid:", ellps);
138 }
139 */
140
141 if (!(print_flag & PRINT_SH)) {
142 fprintf(stdout, "%-*s %s\n", width, "datum:", datum);
143 fprintf(stdout, "%-*s %s\n", width, "ellipsoid:", ellps);
144 }
145 }
146
147 if (print_flag & PRINT_SH) {
148 if (flat_flag)
149 sep=" ";
150 fprintf(stdout, "n=%s%s", north, sep);
151 fprintf(stdout, "s=%s%s", south, sep);
152 fprintf(stdout, "w=%s%s", west, sep);
153 fprintf(stdout, "e=%s%s", east, sep);
154 if (print_flag & PRINT_3D) {
155 fprintf(stdout, "t=%g%s", window->top, sep);
156 fprintf(stdout, "b=%g%s", window->bottom, sep);
157 }
158 fprintf(stdout, "nsres=%s%s", nsres, sep);
159 if (print_flag & PRINT_3D) {
160 fprintf(stdout, "nsres3=%s%s", nsres3, sep);
161 }
162 fprintf(stdout, "ewres=%s%s", ewres, sep);
163 if (print_flag & PRINT_3D) {
164 fprintf(stdout, "ewres3=%s%s", ewres3, sep);
165 fprintf(stdout, "tbres=%s%s", tbres, sep);
166 }
167 fprintf(stdout, "rows=%d%s", window->rows, sep);
168 if (print_flag & PRINT_3D)
169 fprintf(stdout, "rows3=%d%s", window->rows3, sep);
170 fprintf(stdout, "cols=%d%s", window->cols, sep);
171 if (print_flag & PRINT_3D) {
172 fprintf(stdout, "cols3=%d%s", window->cols3, sep);
173 fprintf(stdout, "depths=%d%s", window->depths, sep);
174 }
175#ifdef HAVE_LONG_LONG_INT
176 fprintf(stdout, "cells=%lld%s",
177 (long long)window->rows * window->cols, sep);
178 if (print_flag & PRINT_3D)
179 fprintf(stdout, "cells3=%lld%s",
180 (long long)window->rows3 * window->cols3 *
181 window->depths, sep);
182#else
183 fprintf(stdout, "cells=%ld%s", (long)window->rows * window->cols, sep);
184 if (print_flag & PRINT_3D)
185 fprintf(stdout, "cells3=%ld%s",
186 (long)window->rows3 * window->cols3 * window->depths, sep);
187#endif
188 if (flat_flag)
189 fprintf(stdout, "\n");
190 }
191 else {
192 fprintf(stdout, "%-*s %s\n", width, "north:", north);
193 fprintf(stdout, "%-*s %s\n", width, "south:", south);
194 fprintf(stdout, "%-*s %s\n", width, "west:", west);
195 fprintf(stdout, "%-*s %s\n", width, "east:", east);
196 if (print_flag & PRINT_3D) {
197 fprintf(stdout, "%-*s %.8f\n", width, "top:", window->top);
198 fprintf(stdout, "%-*s %.8f\n", width, "bottom:",
199 window->bottom);
200 }
201 fprintf(stdout, "%-*s %s\n", width, "nsres:", nsres);
202 if (print_flag & PRINT_3D) {
203 fprintf(stdout, "%-*s %s\n", width, "nsres3:", nsres3);
204 }
205 fprintf(stdout, "%-*s %s\n", width, "ewres:", ewres);
206 if (print_flag & PRINT_3D) {
207 fprintf(stdout, "%-*s %s\n", width, "ewres3:", ewres3);
208 fprintf(stdout, "%-*s %s\n", width, "tbres:", tbres);
209 }
210
211 fprintf(stdout, "%-*s %d\n", width, "rows:", window->rows);
212 if (print_flag & PRINT_3D) {
213 fprintf(stdout, "%-*s %d\n", width, "rows3:", window->rows3);
214 }
215 fprintf(stdout, "%-*s %d\n", width, "cols:", window->cols);
216 if (print_flag & PRINT_3D) {
217 fprintf(stdout, "%-*s %d\n", width, "cols3:", window->cols3);
218 fprintf(stdout, "%-*s %d\n", width, "depths:",
219 window->depths);
220 }
221#ifdef HAVE_LONG_LONG_INT
222 fprintf(stdout, "%-*s %lld\n", width, "cells:",
223 (long long)window->rows * window->cols);
224 if (print_flag & PRINT_3D)
225 fprintf(stdout, "%-*s %lld\n", width, "cells3:",
226 (long long)window->rows3 * window->cols3 *
227 window->depths);
228#else
229 fprintf(stdout, "%-*s %ld\n", width, "cells:",
230 (long)window->rows * window->cols);
231 if (print_flag & PRINT_3D)
232 fprintf(stdout, "%-*s %ld\n", width, "cells3:",
233 (long)window->rows3 * window->cols3 * window->depths);
234#endif
235 }
236 }
237
238 /* flag.lprint: show corners and center in lat/long MN 2001 */
239 if (print_flag & PRINT_LL) {
240 double lo1, la1, lo2, la2, lo3, la3, lo4, la4, loc, lac;
241
242 /* if coordinates are not in lat/long format, transform them: */
243 if ((G_projection() != PROJECTION_LL) && window->proj != 0) {
244 /* projection information of input map */
245 struct Key_Value *in_proj_info, *in_unit_info;
246 struct pj_info iproj, oproj, tproj; /* proj parameters */
247
248 /* read current projection info */
249 if ((in_proj_info = G_get_projinfo()) == NULL)
250 G_fatal_error(_("Can't get projection info of current location"));
251
252 if ((in_unit_info = G_get_projunits()) == NULL)
253 G_fatal_error(_("Can't get projection units of current location"));
254
255 if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
256 G_fatal_error(_("Can't get projection key values of current location"));
257
258 G_free_key_value(in_proj_info);
259 G_free_key_value(in_unit_info);
260
261 oproj.pj = NULL;
262 tproj.def = NULL;
263
264 if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
265 G_fatal_error(_("Unable to initialize coordinate transformation"));
266
267 /*
268 * 1 ------ 2
269 * | | map corners
270 * | |
271 * 4--------3
272 */
273
274 latitude = window->north;
275 longitude = window->west;
276 /* get lat/long w/ same datum/ellipsoid as input */
277 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
278 &longitude, &latitude, NULL) < 0)
279 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
280 "GPJ_transform()");
281
282 lo1 = longitude;
283 la1 = latitude;
284
285 latitude = window->north;
286 longitude = window->east;
287 /* get lat/long w/ same datum/ellipsoid as input */
288 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
289 &longitude, &latitude, NULL) < 0)
290 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
291 "GPJ_transform()");
292
293 lo2 = longitude;
294 la2 = latitude;
295
296 latitude = window->south;
297 longitude = window->east;
298 /* get lat/long w/ same datum/ellipsoid as input */
299 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
300 &longitude, &latitude, NULL) < 0)
301 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
302 "GPJ_transform()");
303
304 lo3 = longitude;
305 la3 = latitude;
306
307 latitude = window->south;
308 longitude = window->west;
309 /* get lat/long w/ same datum/ellipsoid as input */
310 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
311 &longitude, &latitude, NULL) < 0)
312 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
313 "GPJ_transform()");
314
315 lo4 = longitude;
316 la4 = latitude;
317
318 /* center coordinates of the current region,
319 * not average of the projected corner coordinates */
320 latitude = (window->north + window->south) / 2.;
321 longitude = (window->west + window->east) / 2.;
322 /* get lat/long w/ same datum/ellipsoid as input */
323 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
324 &longitude, &latitude, NULL) < 0)
325 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
326 "GPJ_transform()");
327
328 loc = longitude;
329 lac = latitude;
330
331 if (print_flag & PRINT_SH) {
332 fprintf(stdout, "nw_long=%.8f\nnw_lat=%.8f\n", lo1, la1);
333 fprintf(stdout, "ne_long=%.8f\nne_lat=%.8f\n", lo2, la2);
334 fprintf(stdout, "se_long=%.8f\nse_lat=%.8f\n", lo3, la3);
335 fprintf(stdout, "sw_long=%.8f\nsw_lat=%.8f\n", lo4, la4);
336 fprintf(stdout, "center_long=%.8f\n", loc);
337 fprintf(stdout, "center_lat=%.8f\n", lac);
338
339 }
340 else {
341 G_format_easting(lo1, buf, PROJECTION_LL);
342 fprintf(stdout, "%-*s long: %s ", width, "north-west corner:",
343 buf);
344 G_format_northing(la1, buf, PROJECTION_LL);
345 fprintf(stdout, "lat: %s\n", buf);
346
347 G_format_easting(lo2, buf, PROJECTION_LL);
348 fprintf(stdout, "%-*s long: %s ", width, "north-east corner:",
349 buf);
350 G_format_northing(la2, buf, PROJECTION_LL);
351 fprintf(stdout, "lat: %s\n", buf);
352
353 G_format_easting(lo3, buf, PROJECTION_LL);
354 fprintf(stdout, "%-*s long: %s ", width, "south-east corner:",
355 buf);
356 G_format_northing(la3, buf, PROJECTION_LL);
357 fprintf(stdout, "lat: %s\n", buf);
358
359 G_format_easting(lo4, buf, PROJECTION_LL);
360 fprintf(stdout, "%-*s long: %s ", width, "south-west corner:",
361 buf);
362 G_format_northing(la4, buf, PROJECTION_LL);
363 fprintf(stdout, "lat: %s\n", buf);
364
365 G_format_easting(loc, buf, PROJECTION_LL);
366 fprintf(stdout, "%-*s %11s\n", width, "center longitude:",
367 buf);
368
369 G_format_northing(lac, buf, PROJECTION_LL);
370 fprintf(stdout, "%-*s %11s\n", width, "center latitude:",
371 buf);
372 }
373
374 if (!(print_flag & PRINT_REG)) {
375 if (print_flag & PRINT_SH) {
376 fprintf(stdout, "rows=%d\n", window->rows);
377 fprintf(stdout, "cols=%d\n", window->cols);
378 }
379 else {
380 fprintf(stdout, "%-*s %d\n", width, "rows:",
381 window->rows);
382 fprintf(stdout, "%-*s %d\n", width, "cols:",
383 window->cols);
384 }
385 }
386 }
387 else { /* in lat/long already */
388
389 if (window->proj != 0)
390 G_message(_("You are already in Lat/Long. Use the -p flag instead."));
391 else
392 G_message(_("You are in a simple XY location, projection to Lat/Lon "
393 "is not possible. Use the -p flag instead."));
394 }
395 }
396
397 /* flag.eprint */
398 if (print_flag & PRINT_EXTENT) {
399 if (print_flag & PRINT_SH) {
400 fprintf(stdout, "ns_extent=%f\n", window->north - window->south);
401 fprintf(stdout, "ew_extent=%f\n", window->east - window->west);
402 }
403 else {
404 if (G_projection() != PROJECTION_LL) {
405 fprintf(stdout, "%-*s %f\n", width, "north-south extent:",
406 window->north - window->south);
407 fprintf(stdout, "%-*s %f\n", width, "east-west extent:",
408 window->east - window->west);
409 }
410 else {
411 G_format_northing(window->north - window->south, buf,
412 PROJECTION_LL);
413 fprintf(stdout, "%-*s %s\n", width, "north-south extent:",
414 buf);
415 G_format_easting(window->east - window->west, buf,
416 PROJECTION_LL);
417 fprintf(stdout, "%-*s %s\n", width, "east-west extent:", buf);
418 }
419 }
420 }
421
422 /* flag.center */
423 if (print_flag & PRINT_CENTER) {
424 if (print_flag & PRINT_SH) {
425 fprintf(stdout, "center_easting=%f\n",
426 (window->west + window->east) / 2.);
427 fprintf(stdout, "center_northing=%f\n",
428 (window->north + window->south) / 2.);
429 }
430 else {
431 if (G_projection() != PROJECTION_LL) {
432 fprintf(stdout, "%-*s %f\n", width, "center easting:",
433 (window->west + window->east) / 2.);
434 fprintf(stdout, "%-*s %f\n", width, "center northing:",
435 (window->north + window->south) / 2.);
436 }
437 else {
438 G_format_northing((window->north + window->south) / 2.,
439 buf, PROJECTION_LL);
440 fprintf(stdout, "%-*s %s\n", width, "north-south center:",
441 buf);
442 G_format_easting((window->west + window->east) / 2.,
443 buf, PROJECTION_LL);
444 fprintf(stdout, "%-*s %s\n", width, "east-west center:", buf);
445 }
446 }
447 }
448
449
450 /* flag.gmt_style */
451 if (print_flag & PRINT_GMT)
452 fprintf(stdout, "%s/%s/%s/%s\n", west, east, south, north);
453
454 /* flag.wms_style */
455 if (print_flag & PRINT_WMS) {
456 G_format_northing(window->north, north, -1);
457 G_format_northing(window->south, south, -1);
458 G_format_easting(window->east, east, -1);
459 G_format_easting(window->west, west, -1);
460 fprintf(stdout, "bbox=%s,%s,%s,%s\n", west, south, east, north);
461 }
462
463
464 /* flag.nangle */
465 if (print_flag & PRINT_NANGLE) {
466 double convergence;
467
468 if (G_projection() == PROJECTION_XY)
469 convergence = 0./0.;
470 else if (G_projection() == PROJECTION_LL)
471 convergence = 0.0;
472 else {
473 struct Key_Value *in_proj_info, *in_unit_info;
474 /* proj parameters */
475 struct pj_info iproj, oproj, tproj;
476#ifdef HAVE_PROJ_H
477 PJ_COORD c;
478 PJ_FACTORS fact;
479#else
480 struct FACTORS fact;
481 LP lp;
482#endif
483
484 /* read current projection info */
485 if ((in_proj_info = G_get_projinfo()) == NULL)
486 G_fatal_error(_("Can't get projection info of current location"));
487
488 if ((in_unit_info = G_get_projunits()) == NULL)
489 G_fatal_error(_("Can't get projection units of current location"));
490
491 if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
492 G_fatal_error(_("Can't get projection key values of current location"));
493
494 G_free_key_value(in_proj_info);
495 G_free_key_value(in_unit_info);
496
497 oproj.pj = NULL;
498 tproj.def = NULL;
499
500 if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
501 G_fatal_error(_("Unable to initialize coordinate transformation"));
502
503 /* center coordinates of the current region,
504 * not average of the projected corner coordinates */
505 latitude = (window->north + window->south) / 2.;
506 longitude = (window->west + window->east) / 2.;
507 /* get lat/long w/ same datum/ellipsoid as input */
508 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
509 &longitude, &latitude, NULL) < 0)
510 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
511 "GPJ_transform()");
512
513#ifdef HAVE_PROJ_H
514 c.lpzt.lam = DEG2RAD(longitude);
515 c.lpzt.phi = DEG2RAD(latitude);
516 c.lpzt.z = 0;
517 c.lpzt.t = 0;
518 fact = proj_factors(iproj.pj, c);
519 convergence = RAD2DEG(fact.meridian_convergence);
520#else
521 G_zero(&fact, sizeof(struct FACTORS));
522 lp.u = DEG2RAD(longitude);
523 lp.v = DEG2RAD(latitude);
524 pj_factors(lp, iproj.pj, 0.0, &fact);
525 convergence = RAD2DEG(fact.conv);
526#endif
527 }
528
529 if (print_flag & PRINT_SH)
530 fprintf(stdout, "converge_angle=%f\n", convergence);
531 else
532 fprintf(stdout, "%-*s %f\n", width, "convergence angle:",
533 convergence);
534 }
535
536 /* flag.bbox
537 Calculate the largest bounding box in lat-lon coordinates
538 and print it to stdout
539 */
540 if (print_flag & PRINT_MBBOX) {
541 double sh_ll_w, sh_ll_e, sh_ll_n, sh_ll_s, loc;
542
543 /* Needed to calculate the LL bounding box */
544 if ((G_projection() != PROJECTION_XY)) {
545 /* projection information of input and output map */
546 struct Key_Value *in_proj_info, *in_unit_info, *out_proj_info,
547 *out_unit_info;
548 struct pj_info iproj; /* input map proj parameters */
549 struct pj_info oproj; /* output map proj parameters */
550 struct pj_info tproj; /* transformation parameters */
551 char buff[100], dum[100];
552 int r, c;
553
554 /* read current projection info */
555 if ((in_proj_info = G_get_projinfo()) == NULL)
556 G_fatal_error(_("Can't get projection info of current location"));
557 /* do not wrap to -180, 180, otherwise east can be < west */
558 G_set_key_value("+over", "defined", in_proj_info);
559
560 if ((in_unit_info = G_get_projunits()) == NULL)
561 G_fatal_error(_("Can't get projection units of current location"));
562
563 if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
564 G_fatal_error(_("Can't get projection key values of current location"));
565
566 /* output projection to lat/long and wgs84 ellipsoid */
567 out_proj_info = G_create_key_value();
568 out_unit_info = G_create_key_value();
569
570 G_set_key_value("proj", "ll", out_proj_info);
571
572 if (G_get_datumparams_from_projinfo(in_proj_info, buff, dum) < 0)
573 G_fatal_error(_("WGS84 output not possible as this location does not contain "
574 "datum transformation parameters. Try running g.setproj."));
575 else
576 G_set_key_value("datum", "wgs84", out_proj_info);
577
578 G_set_key_value("unit", "degree", out_unit_info);
579 G_set_key_value("units", "degrees", out_unit_info);
580 G_set_key_value("meters", "1.0", out_unit_info);
581
582 if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
583 G_fatal_error(_("Unable to update lat/long projection parameters"));
584
585 G_free_key_value(in_proj_info);
586 G_free_key_value(in_unit_info);
587 G_free_key_value(out_proj_info);
588 G_free_key_value(out_unit_info);
589
590 tproj.def = NULL;
591
592 if (GPJ_init_transform(&iproj, &oproj, &tproj) < 0)
593 G_fatal_error(_("Unable to initialize coordinate transformation"));
594
595 /*Calculate the largest bounding box */
596
597 /* center */
598 latitude = (window->north + window->south) / 2.;
599 longitude = (window->west + window->east) / 2.;
600 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
601 &longitude, &latitude, NULL) < 0)
602 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
603 "GPJ_transform()");
604
605 sh_ll_w = sh_ll_e = longitude;
606 sh_ll_n = sh_ll_s = latitude;
607
608 /* western and eastern border */
609 for (r = 0; r <= window->rows; r++) {
610 latitude = window->north - r * window->ns_res;
611 if (r == window->rows)
612 latitude = window->south;
613 longitude = window->west;
614 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
615 &longitude, &latitude, NULL) < 0)
616 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
617 "GPJ_transform()");
618
619 if (sh_ll_n < latitude)
620 sh_ll_n = latitude;
621 if (sh_ll_s > latitude)
622 sh_ll_s = latitude;
623
624 if (sh_ll_e < longitude)
625 sh_ll_e = longitude;
626 if (sh_ll_w > longitude)
627 sh_ll_w = longitude;
628
629 latitude = window->north - r * window->ns_res;
630 if (r == window->rows)
631 latitude = window->south;
632 longitude = window->east;
633 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
634 &longitude, &latitude, NULL) < 0)
635 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
636 "GPJ_transform()");
637
638 if (sh_ll_n < latitude)
639 sh_ll_n = latitude;
640 if (sh_ll_s > latitude)
641 sh_ll_s = latitude;
642
643 if (sh_ll_e < longitude)
644 sh_ll_e = longitude;
645 if (sh_ll_w > longitude)
646 sh_ll_w = longitude;
647 }
648
649 /* northern and southern border */
650 for (c = 1; c < window->cols; c++) {
651 latitude = window->north;
652 longitude = window->west + c * window->ew_res;
653 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
654 &longitude, &latitude, NULL) < 0)
655 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
656 "GPJ_transform()");
657
658 if (sh_ll_n < latitude)
659 sh_ll_n = latitude;
660 if (sh_ll_s > latitude)
661 sh_ll_s = latitude;
662
663 if (sh_ll_e < longitude)
664 sh_ll_e = longitude;
665 if (sh_ll_w > longitude)
666 sh_ll_w = longitude;
667
668 latitude = window->south;
669 longitude = window->west + c * window->ew_res;
670 if (GPJ_transform(&iproj, &oproj, &tproj, PJ_FWD,
671 &longitude, &latitude, NULL) < 0)
672 G_fatal_error(_("Error in %s (projection of input coordinate pair)"),
673 "GPJ_transform()");
674
675 if (sh_ll_n < latitude)
676 sh_ll_n = latitude;
677 if (sh_ll_s > latitude)
678 sh_ll_s = latitude;
679
680 if (sh_ll_e < longitude)
681 sh_ll_e = longitude;
682 if (sh_ll_w > longitude)
683 sh_ll_w = longitude;
684 }
685
686 loc = (sh_ll_e + sh_ll_w) / 2.;
687 loc += get_shift(loc);
688 sh_ll_w += get_shift(sh_ll_w);
689 sh_ll_e += get_shift(sh_ll_e);
690
691 /* print the largest bounding box */
692 if (print_flag & PRINT_SH) {
693 fprintf(stdout, "ll_n=%.8f\n", sh_ll_n);
694 fprintf(stdout, "ll_s=%.8f\n", sh_ll_s);
695 fprintf(stdout, "ll_w=%.8f\n", sh_ll_w);
696 fprintf(stdout, "ll_e=%.8f\n", sh_ll_e);
697 /* center of the largest bounding box */
698 fprintf(stdout, "ll_clon=%.8f\n", loc);
699 fprintf(stdout, "ll_clat=%.8f\n",
700 (sh_ll_n + sh_ll_s) / 2.);
701 }
702 else {
703 G_format_northing(sh_ll_n, buf, PROJECTION_LL);
704 fprintf(stdout, "%-*s %s\n", width, "north latitude:", buf);
705 G_format_northing(sh_ll_s, buf, PROJECTION_LL);
706 fprintf(stdout, "%-*s %s\n", width, "south latitude:", buf);
707 G_format_easting(sh_ll_w, buf, PROJECTION_LL);
708 fprintf(stdout, "%-*s %s\n", width, "west longitude:", buf);
709 G_format_easting(sh_ll_e, buf, PROJECTION_LL);
710 fprintf(stdout, "%-*s %s\n", width, "east longitude:", buf);
711 /* center of the largest bounding box */
712 G_format_easting(loc, buf, PROJECTION_LL);
713 fprintf(stdout, "%-*s %s\n", width, "center longitude:", buf);
714 G_format_northing((sh_ll_n + sh_ll_s) / 2., buf,
715 PROJECTION_LL);
716 fprintf(stdout, "%-*s %s\n", width, "center latitude:", buf);
717 }
718
719 /*It should be calculated which number of rows and cols we have in LL */
720 /*
721 fprintf (stdout, "LL_ROWS=%f \n",sh_ll_rows);
722 fprintf (stdout, "LL_COLS=%f \n",sh_ll_cols);
723 */
724 }
725 else {
726 G_warning(_("Lat/Long calculations are not possible from a simple XY system"));
727 }
728 }
729}
Note: See TracBrowser for help on using the repository browser.