[GRASS-SVN] r43941 -
grass/branches/releasebranch_6_4/raster/r.proj.seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Oct 17 03:03:35 EDT 2010
Author: hamish
Date: 2010-10-17 00:03:35 -0700 (Sun, 17 Oct 2010)
New Revision: 43941
Modified:
grass/branches/releasebranch_6_4/raster/r.proj.seg/description.html
grass/branches/releasebranch_6_4/raster/r.proj.seg/main.c
Log:
backport -p,-g print source region flags from devbr6
Modified: grass/branches/releasebranch_6_4/raster/r.proj.seg/description.html
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.proj.seg/description.html 2010-10-17 07:01:41 UTC (rev 43940)
+++ grass/branches/releasebranch_6_4/raster/r.proj.seg/description.html 2010-10-17 07:03:35 UTC (rev 43941)
@@ -137,8 +137,19 @@
To avoid excessive time consumption when reprojecting a map the region and
resolution of the target location should be set appropriately beforehand.
-A simple way to do this is to generate a vector "box" map of the region in
-the source location using <em><a href="v.in.region.html">v.in.region</a></em>.
+<p>
+A simple way to do this is to check the projected bounds of the input map
+in the current location's projection using the <b>-p</b> flag. The <b>-g</b>
+flag reports the same thing, but in a form which can be directly cut and
+pasted into a <em>g.region</em> command. After setting the region in that
+way you might check the cell resolution with "<em>g.region -p</em>" then
+snap it to a regular grid with <em>g.region</em>'s -a flag. E.g.
+<tt>g.region -a res=5 -p</tt>. Note that this is just a rough guide.
+
+<p>
+A more involved, but more accurate, way to do this is to generate a vector
+"box" map of the region in the source location using
+ <em><a href="v.in.region.html">v.in.region</a></em>.
This "box" map is then reprojected into the target location with
<em><a href="v.proj.html">v.proj</a></em>.
Next the region in the target location is set to the extent of the new vector
@@ -157,9 +168,66 @@
<h2>EXAMPLES</h2>
-<h3>v.in.region method</h3>
+<h4>Inline method</h4>
+With GRASS running in the destination location use the <b>-g</b> flag to
+show the input map's bounds once reprojected into the current map projection,
+then use that to set the region bounds before performing the reprojection:
+
<div class="code"><pre>
+# calculate where output map will be
+GRASS> r.proj input=elevation location=ll_wgs84 mapset=user1 -p
+Source cols: 8162
+Source rows: 12277
+Local north: -4265502.30382993
+Local south: -4473453.15255565
+Local west: 14271663.19157564
+Local east: 14409956.2693866
+# same calculation, but in a form which can be cut and pasted into a g.region call
+GRASS> r.proj input=elevation location=ll_wgs84 mapset=user1 -g
+n=-4265502.30382993 s=-4473453.15255565 w=14271663.19157564 e=14409956.2693866 rows=12277 col
+s=8162
+
+GRASS> g.region n=-4265502.30382993 s=-4473453.15255565 \
+ w=14271663.19157564 e=14409956.2693866 rows=12277 cols=8162 -p
+projection: 99 (Mercator)
+zone: 0
+datum: wgs84
+ellipsoid: wgs84
+north: -4265502.30382993
+south: -4473453.15255565
+west: 14271663.19157564
+east: 14409956.2693866
+nsres: 16.93824621
+ewres: 16.94352828
+rows: 12277
+cols: 8162
+cells: 100204874
+
+# round resolution to something cleaner
+GRASS> g.region res=17 -a -p
+projection: 99 (Mercator)
+zone: 0
+datum: wgs84
+ellipsoid: wgs84
+north: -4265487
+south: -4473465
+west: 14271653
+east: 14409965
+nsres: 17
+ewres: 17
+rows: 12234
+cols: 8136
+cells: 99535824
+
+# finally, perform the reprojection
+GRASS> r.proj input=elevation location=ll_wgs84 mapset=user1 memory=800
+</pre></div>
+
+
+<h4>v.in.region method</h4>
+<div class="code"><pre>
+
# In the source location, use v.in.region to generate a bounding box around the
# region of interest:
Modified: grass/branches/releasebranch_6_4/raster/r.proj.seg/main.c
===================================================================
--- grass/branches/releasebranch_6_4/raster/r.proj.seg/main.c 2010-10-17 07:01:41 UTC (rev 43940)
+++ grass/branches/releasebranch_6_4/raster/r.proj.seg/main.c 2010-10-17 07:03:35 UTC (rev 43941)
@@ -87,7 +87,8 @@
row, col, /* counters */
irows, icols, /* original rows, cols */
orows, ocols, have_colors, /* Input map has a colour table */
- overwrite; /* Overwrite */
+ overwrite, /* Overwrite */
+ curr_proj; /* output projection (see gis.h) */
void *obuffer, /* buffer that holds one output row */
*obufptr; /* column ptr in output buffer */
@@ -100,6 +101,7 @@
row_idx, /* row index in input matrix */
onorth, osouth, /* save original border coords */
oeast, owest, inorth, isouth, ieast, iwest;
+ char north_str[30], south_str[30], east_str[30], west_str[30];
struct Colors colr; /* Input map colour table */
struct History history;
@@ -114,7 +116,9 @@
struct GModule *module;
struct Flag *list, /* list files in source location */
- *nocrop; /* don't crop output map */
+ *nocrop, /* don't crop output map */
+ *print_bounds, /* print output bounds and exit */
+ *gprint_bounds; /* same but print shell style */
struct Option *imapset, /* name of input mapset */
*inmap, /* name of input layer */
@@ -192,6 +196,17 @@
nocrop->key = 'n';
nocrop->description = _("Do not perform region cropping optimization");
+ print_bounds = G_define_flag();
+ print_bounds->key = 'p';
+ print_bounds->description =
+ _("Print input map's bounds in the current projection and exit");
+
+ gprint_bounds = G_define_flag();
+ gprint_bounds->key = 'g';
+ gprint_bounds->description =
+ _("Print input map's bounds in the current projection and exit (shell style)");
+
+
/* The parser checks if the map already exists in current mapset,
we switch out the check and do it
in the module after the parser */
@@ -200,6 +215,7 @@
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+
/* get the method */
for (method = 0; (ipolname = menu[method].name); method++)
if (strcmp(ipolname, interpol->answer) == 0)
@@ -222,6 +238,10 @@
G_get_window(&outcellhd);
+ if(gprint_bounds->answer && !print_bounds->answer)
+ print_bounds->answer = gprint_bounds->answer;
+ curr_proj = G_projection();
+
/* Get projection info for output mapset */
if ((out_proj_info = G_get_projinfo()) == NULL)
G_fatal_error(_("Unable to get projection info of output raster map"));
@@ -302,6 +322,40 @@
orows = outcellhd.rows;
ocols = outcellhd.cols;
+
+ if (print_bounds->answer) {
+ G_message(_("Input map <%s@%s> in location <%s>:"),
+ inmap->answer, setname, inlocation->answer);
+
+ if (pj_do_proj(&iwest, &isouth, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+ if (pj_do_proj(&ieast, &inorth, &iproj, &oproj) < 0)
+ G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+
+ G_format_northing(inorth, north_str, curr_proj);
+ G_format_northing(isouth, south_str, curr_proj);
+ G_format_easting(ieast, east_str, curr_proj);
+ G_format_easting(iwest, west_str, curr_proj);
+
+ if(gprint_bounds->answer) {
+ fprintf(stdout, "n=%s s=%s w=%s e=%s rows=%d cols=%d\n",
+ north_str, south_str, west_str, east_str, irows, icols);
+ }
+ else {
+ fprintf(stdout, "Source cols: %d\n", icols);
+ fprintf(stdout, "Source rows: %d\n", irows);
+ fprintf(stdout, "Local north: %s\n", north_str);
+ fprintf(stdout, "Local south: %s\n", south_str);
+ fprintf(stdout, "Local west: %s\n", west_str);
+ fprintf(stdout, "Local east: %s\n", east_str);
+ }
+
+ /* somehow approximate local ewres, nsres ?? (use 'g.region -m' on lat/lon side) */
+
+ exit(EXIT_SUCCESS);
+ }
+
+
/* Cut non-overlapping parts of input map */
if (!nocrop->answer)
bordwalk(&outcellhd, &incellhd, &oproj, &iproj);
More information about the grass-commit
mailing list