[GRASS-SVN] r33995 - grass/branches/develbranch_6/general/g.region
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Oct 23 21:18:34 EDT 2008
Author: hamish
Date: 2008-10-23 21:18:33 -0400 (Thu, 23 Oct 2008)
New Revision: 33995
Modified:
grass/branches/develbranch_6/general/g.region/local_proto.h
grass/branches/develbranch_6/general/g.region/main.c
grass/branches/develbranch_6/general/g.region/printwindow.c
Log:
add -n flag to output convergence angle
Modified: grass/branches/develbranch_6/general/g.region/local_proto.h
===================================================================
--- grass/branches/develbranch_6/general/g.region/local_proto.h 2008-10-23 21:29:19 UTC (rev 33994)
+++ grass/branches/develbranch_6/general/g.region/local_proto.h 2008-10-24 01:18:33 UTC (rev 33995)
@@ -9,6 +9,7 @@
#define PRINT_METERS 0x20
#define PRINT_3D 0x40
#define PRINT_MBBOX 0x80
+#define PRINT_NANGLE 0x100
/* adjust.c */
int adjust_window(struct Cell_head *, int, int, int);
Modified: grass/branches/develbranch_6/general/g.region/main.c
===================================================================
--- grass/branches/develbranch_6/general/g.region/main.c 2008-10-23 21:29:19 UTC (rev 33994)
+++ grass/branches/develbranch_6/general/g.region/main.c 2008-10-24 01:18:33 UTC (rev 33995)
@@ -45,11 +45,7 @@
struct
{
struct Flag
- *update,
- *print,
- *gprint,
- *lprint,
- *eprint,
+ *update, *print, *gprint, *lprint, *eprint, *nangle,
*center, *res_set, *dist_res, *dflt, *z, *savedefault, *bbox;
} flag;
struct
@@ -113,6 +109,14 @@
_("Print region resolution in meters (geodesic)");
flag.dist_res->guisection = _("Print");
+ flag.nangle = G_define_flag();
+ flag.nangle->key = 'n';
+ flag.nangle->label = _("Print the convergence angle (degrees CCW)");
+ flag.nangle->description =
+ _("The difference between the projection's grid north and true north, "
+ "measured at the center coordinates of the current region.");
+ flag.nangle->guisection = _("Print");
+
flag.z = G_define_flag();
flag.z->key = '3';
flag.z->description = _("Print also 3D settings");
@@ -366,6 +370,9 @@
if (flag.center->answer)
print_flag |= PRINT_CENTER;
+ if (flag.nangle->answer)
+ print_flag |= PRINT_NANGLE;
+
if (flag.dist_res->answer)
print_flag |= PRINT_METERS;
Modified: grass/branches/develbranch_6/general/g.region/printwindow.c
===================================================================
--- grass/branches/develbranch_6/general/g.region/printwindow.c 2008-10-23 21:29:19 UTC (rev 33994)
+++ grass/branches/develbranch_6/general/g.region/printwindow.c 2008-10-24 01:18:33 UTC (rev 33995)
@@ -1,10 +1,15 @@
#include <string.h>
#include <stdlib.h>
+#include <projects.h>
#include <grass/gis.h>
#include <grass/gprojects.h>
#include <grass/glocale.h>
#include "local_proto.h"
+#define DEG2RAD(a) ((a) * M_PI / 180.0)
+#define RAD2DEG(a) ((a) * 180.0 / M_PI)
+
+
int print_window(struct Cell_head *window, int print_flag)
{
char *prj, *datum, *ellps;
@@ -53,7 +58,7 @@
if (print_flag & PRINT_CENTER || print_flag & PRINT_MBBOX)
width = 16;
- if (print_flag & PRINT_LL)
+ if (print_flag & PRINT_LL || print_flag & PRINT_NANGLE)
width = 18;
if (print_flag & PRINT_EXTENT)
@@ -463,6 +468,122 @@
}
}
+ /* flag.nangle */
+ if (print_flag & PRINT_NANGLE) {
+ double convergence;
+
+ if (G_projection() == PROJECTION_XY)
+ convergence = 0./0.;
+ else if (G_projection() == PROJECTION_LL)
+ convergence = 0.0;
+ else {
+ struct Key_Value *in_proj_info, *in_unit_info;
+ double lo1, la1, lo2, la2, lo3, la3, lo4, la4;
+ double mid_n_lo, mid_n_la, mid_s_lo, mid_s_la;
+ double lat_center, lon_center;
+ struct pj_info iproj, oproj; /* proj parameters */
+ struct FACTORS fact;
+
+ /* read current projection info */
+ if ((in_proj_info = G_get_projinfo()) == NULL)
+ G_fatal_error(_("Can't get projection info of current location"));
+
+ if ((in_unit_info = G_get_projunits()) == NULL)
+ G_fatal_error(_("Can't get projection units of current location"));
+
+ if (pj_get_kv(&iproj, in_proj_info, in_unit_info) < 0)
+ G_fatal_error(_("Can't get projection key values of current location"));
+
+ /* output projection to lat/long w/ same ellipsoid as input */
+ oproj.zone = 0;
+ oproj.meters = 1.;
+ sprintf(oproj.proj, "ll");
+ if ((oproj.pj = pj_latlong_from_proj(iproj.pj)) == NULL)
+ G_fatal_error(_("Unable to update lat/long projection parameters"));
+
+ /* for DEBUG
+ pj_print_proj_params(&iproj, &oproj);
+ */
+
+ /* do the transform
+ * syntax: pj_do_proj(outx, outy, in_info, out_info)
+ *
+ * 1 ------ 2
+ * | | map corners
+ * | |
+ * 4--------3
+ */
+
+ latitude = window->north;
+ longitude = window->west;
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+
+ lo1 = longitude;
+ la1 = latitude;
+
+ latitude = window->north;
+ longitude = window->east;
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+
+ lo2 = longitude;
+ la2 = latitude;
+
+ latitude = window->south;
+ longitude = window->east;
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+
+ lo3 = longitude;
+ la3 = latitude;
+
+ latitude = window->south;
+ longitude = window->west;
+ if (pj_do_proj(&longitude, &latitude, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+
+ lo4 = longitude;
+ la4 = latitude;
+
+ /*
+ * map corners and side mids:
+ * mid_n
+ * 1 ---+---2
+ * | |
+ * mid_w + + mid_e
+ * | |
+ * 4----+---3
+ * mid_s
+ *
+ * lo: longitude
+ * la: latitude
+ */
+
+ /* side mids for easting, northing center: */
+ mid_n_lo = (lo2 + lo1) / 2.;
+ mid_n_la = (la2 + la1) / 2.;
+ mid_s_lo = (lo3 + lo4) / 2.;
+ mid_s_la = (la3 + la4) / 2.;
+
+ lat_center = (mid_n_la + mid_s_la) / 2.;
+ lon_center = (mid_n_lo + mid_s_lo) / 2.;
+
+ LP lp = { DEG2RAD(lon_center), DEG2RAD(lat_center) };
+ pj_factors(lp, iproj.pj, 0.0, &fact);
+ convergence = RAD2DEG(fact.conv);
+
+ G_free_key_value(in_proj_info);
+ G_free_key_value(in_unit_info);
+ }
+
+ if (print_flag & PRINT_SH)
+ fprintf(stdout, "converge_angle=%f\n", convergence);
+ else
+ fprintf(stdout, "%-*s %f\n", width, "convergence angle:",
+ convergence);
+ }
+
/* flag.bbox
Calculate the largest boudingbox in lat-lon coordinates
and print it to stdout
More information about the grass-commit
mailing list