[GRASS-SVN] r30699 - grass/trunk/general/g.region
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Mar 22 16:39:13 EDT 2008
Author: martinl
Date: 2008-03-22 16:39:13 -0400 (Sat, 22 Mar 2008)
New Revision: 30699
Added:
grass/trunk/general/g.region/adjust.c
grass/trunk/general/g.region/description.html
grass/trunk/general/g.region/local_proto.h
grass/trunk/general/g.region/main.c
grass/trunk/general/g.region/printwindow.c
grass/trunk/general/g.region/zoom.c
Removed:
grass/trunk/general/g.region/cmd/
Modified:
grass/trunk/general/g.region/Makefile
Log:
g.region: source code moved from 'cmd' directory to . ('cmd' directory not needed here)
Modified: grass/trunk/general/g.region/Makefile
===================================================================
--- grass/trunk/general/g.region/Makefile 2008-03-22 20:32:39 UTC (rev 30698)
+++ grass/trunk/general/g.region/Makefile 2008-03-22 20:39:13 UTC (rev 30699)
@@ -1,12 +1,16 @@
-
MODULE_TOPDIR = ../..
-SUBDIRS = \
- cmd
+PGM = g.region
-include $(MODULE_TOPDIR)/include/Make/Dir.make
+LIBES = $(VECTLIB) $(G3DLIB) $(GPROJLIB) $(GISLIB) $(DATETIMELIB)
+DEPENDENCIES = $(VECTDEP) $(G3DDEP) $(GPROJDEP) $(GISDEP) $(DATETIMEDEP)
+EXTRA_INC = $(VECT_INC) $(PROJINC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
-default: parsubdirs
+include $(MODULE_TOPDIR)/include/Make/Module.make
-clean: cleansubdirs
+default: cmd
+
+
+
Copied: grass/trunk/general/g.region/adjust.c (from rev 30698, grass/trunk/general/g.region/cmd/adjust.c)
===================================================================
--- grass/trunk/general/g.region/adjust.c (rev 0)
+++ grass/trunk/general/g.region/adjust.c 2008-03-22 20:39:13 UTC (rev 30699)
@@ -0,0 +1,15 @@
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+int
+adjust_window (struct Cell_head *window,int row_flag,int col_flag,int depth_flag)
+{
+ char *err = G_adjust_Cell_head3(window, row_flag, col_flag, depth_flag);
+
+ if (err)
+ {
+ G_fatal_error (_("Invalid region: %s"), err);
+ }
+
+ return 0;
+}
Copied: grass/trunk/general/g.region/description.html (from rev 30698, grass/trunk/general/g.region/cmd/description.html)
===================================================================
--- grass/trunk/general/g.region/description.html (rev 0)
+++ grass/trunk/general/g.region/description.html 2008-03-22 20:39:13 UTC (rev 30699)
@@ -0,0 +1,497 @@
+<h2>DESCRIPTION</h2>
+
+The <em>g.region</em> module allows the user to manage the
+settings of the current geographic region. These regional
+boundaries can be set by the user directly and/or set from
+a region definition file (stored under the
+<kbd>windows</kbd> directory in the user's current
+mapset). The user can create, modify, and store as many
+geographic region definitions as desired for any given
+mapset. However, only one of these geographic region
+definitions will be current at any given moment, for a
+specified mapset; i.e., GRASS programs that respect the
+geographic region settings will use the current geographic
+region settings.
+
+
+<h2>DEFINITIONS</h2>
+
+<dl>
+<dt><b>Region:</b>
+
+<dd>In GRASS, a <em>region</em> refers to a geographic area
+with some defined boundaries, based on a specific map
+coordinate system and map projection. Each region also has
+associated with it the specific east-west and north-south
+resolutions of its smallest units (rectangular units called
+"cells").
+
+<p>
+
+The region's boundaries are given as the northernmost,
+southernmost, easternmost, and westernmost points that
+define its extent (cell edges). The north and south boundaries
+are commonly called <em>northings</em>, while the east and west
+boundaries are called <em>eastings</em>.
+
+<p>
+
+The region's cell resolution defines the size of the
+smallest piece of data recognized (imported, analyzed,
+displayed, stored, etc.) by GRASS modules affected by the
+current region settings. The north-south and east-west cell
+resolutions need not be the same, thus allowing non-square
+data cells to exist.
+
+<P>
+Typically all raster and display modules are affected by the current
+region settings, but not vector modules.
+Some special modules diverge from this rule, for example raster import
+modules and <em>v.in.region</em>.
+
+
+<dt><b>Default Region:</b>
+
+<dd>Each GRASS LOCATION has a fixed
+geographic region, called the default geographic region
+(stored in the region file <kbd>DEFAULT_WIND</kbd> under
+the special mapset <kbd>PERMANENT</kbd>), that defines the
+extent of the data base. While this provides a starting
+point for defining new geographic regions, user-defined
+geographic regions need not fall within this geographic
+region. The current region can be reset to the default region
+with the <b>-d</b> flag. The default region is initially set
+when the location is first created and can be reset using the
+<b>-s</b> flag.
+
+<dt><b>Current Region:</b>
+
+<dd>Each mapset has a current geographic region. This
+region defines the geographic area in which all GRASS
+displays and raster analyses will be done. Raster data will be
+resampled, if necessary, to meet the cell resolutions of
+the current geographic region setting.
+
+<dt><b>Saved Regions:</b>
+
+<dd>Each GRASS MAPSET may contain any number of
+pre-defined, and named, geographic regions. These region
+definitions are stored in the user's current mapset
+location under the <kbd>windows</kbd> directory (also
+referred to as the user's saved region definitions).
+Any of these pre-defined geographic regions
+may be selected, by name, to become the current geographic
+region. Users may also access saved region definitions
+stored under other mapsets in the current location, if
+these mapsets are included in the user's mapset search
+path or the '@' operator is used (<tt>region_name at mapset</tt>).
+</dl>
+
+
+<h2>NOTES</h2>
+
+After all updates have been applied, the current region's
+southern and western boundaries are (silently) adjusted so
+that the north/south distance is a multiple of the
+north/south resolution and that the east/west distance is a
+multiple of the east/west resolution.
+
+<p>
+With the <b>-a</b> flag all four boundaries are adjusted
+to be even multiples of the resolution, aligning the region to the
+resolution supplied by the user. The default is to
+align the region resolution to match the region boundaries.
+
+<p>
+The <b>-m</b> flag will report the region resolution in meters. The
+resolution is calculated by averaging the resolution at the region
+boundaries. This resolution is calculated by dividing the geodesic
+distance in meters at the boundary by the number of rows or columns.
+For example the east / west resolution (ewres) is determined from an
+average of the geodesic distances at the North and South boundaries
+divided by the number of columns.
+<!-- add'l info. include?
+Print the region resolution in meters (from geodesic). With no other
+flags the default output format is shell stype (-g). The region resolution
+represents the center of the map. The resolutions are calculated at the four
+outside edges, then the two NS edges are averaged and the two EW edges are
+averaged, the results finally printed.
+-->
+
+
+<p>
+The <b>-p</b> (or <b>-g</b>) option is recognized
+last. This means that all changes are applied to the
+region settings before printing occurs.
+<P>
+The <b>-g</b> flag prints the current region settings in shell script style.
+This format can be given back to <em>g.region</em> on its command line.
+This may also be used to save region settings as shell environment variables
+with the UNIX eval command, "<tt>eval `g.region -g`</tt>".
+
+
+<h3>Additional parameter information:</h3>
+
+<dl>
+<dt><b>3dview=</b><em>name</em>
+<dd>Make current region settings same as those in the named
+3dview file, which holds the region that was current when
+the 3dview was saved.
+
+
+<dt><b>zoom=</b><em>name</em>
+<dd>Shrink current region settings to the smallest region
+encompassing all non-NULL data in the named raster map
+layer that fall inside the user's current region. In this
+way you can tightly zoom in on isolated clumps within a
+bigger map.
+<p>
+If the user also includes the <b>rast=</b><em>name</em>
+option on the command line, <b>zoom=</b><em>name</em> will
+set the current region settings to the smallest region
+encompassing all non-NULL data in the named <b>zoom</b> map
+that fall inside the region stated in the cell header for
+the named <b>raster</b> map.
+
+
+<dt><b>align=</b><em>name</em>
+
+<dd>Set the current resolution equal to that of the named
+raster map, and align the current region to a row and
+column edge in the named map. Alignment only moves the
+existing region edges outward to the edges of the next
+nearest cell in the named raster map - not to the named
+map's edges. To perform the latter function, use the
+<b>rast=</b><em>name</em> option.
+</dl>
+
+
+
+<h2>EXAMPLES</h2>
+
+<dl>
+<dt>
+<div class="code"><pre>
+g.region n=7360100 e=699000
+</pre></div>
+<dd> will reset the northing and easting for the current
+region, but leave the south edge, west edge, and the region
+cell resolutions unchanged.
+
+<p>
+
+<dt><div class="code"><pre>
+g.region n=51:36:05N e=10:10:05E s=51:29:55N w=9:59:55E res=0:00:01
+</pre></div>
+
+<dd> will reset the northing, easting, southing, westing and resolution
+for the current region, here in DMS latitude-longitude style
+(decimal degrees and degrees with decimal minutes can also be used).
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -dp s=698000
+</pre></div>
+
+<dd> will set the current region from the default region
+for the GRASS data base location, reset the south edge to
+698000, and then print the result.
+
+<p>
+
+<dt><div class="code"><pre>
+g.region n=n+1000 w=w-500
+</pre></div>
+
+<dd> The n=<em>value</em> may also be specified as a
+function of its current value: n=n+<em>value</em>
+increases the current northing, while n=n-<em>value</em>
+decreases it. This is also true for s=<em>value</em>,
+e=<em>value</em>, and w=<em>value</em>. In this example
+the current region's northern boundary is extended by 1000
+units and the current region's western boundary is
+decreased by 500 units.
+
+<p>
+
+<dt><div class="code"><pre>
+g.region n=s+1000 e=w+1000
+</pre></div>
+
+<dd> This form allows the user to set the region boundary
+values relative to one another. Here, the northern
+boundary coordinate is set equal to 1000 units larger than
+the southern boundary's coordinate value, and the eastern
+boundary's coordinate value is set equal to 1000 units
+larger than the western boundary's coordinate value. The
+corresponding forms s=n-<em>value</em> and
+
+<p>
+
+w=e-<em>value</em> may be used to set the values of the
+region's southern and western boundaries, relative to the
+northern and eastern boundary values.
+
+
+<dt><div class="code"><pre>
+g.region rast=soils
+</pre></div>
+
+<dd> This form will make the current region settings
+exactly the same as those given in the cell header file for
+the raster map layer <em>soils</em>.
+
+<p>
+
+
+<dt><div class="code"><pre>
+g.region rast=soils zoom=soils
+</pre></div>
+
+<dd> This form will first look up the cell header file for
+the raster map layer <em>soils</em>, use this as the
+current region setting, and then shrink the region down to
+the smallest region which still encompasses all non-NULL
+data in the map layer <em>soils</em>. Note that if the
+parameter <em>rast=soils</em> were not specified, the
+zoom would shrink to encompass all non-NULL data values in
+the soils map that were located within the <i>current region</i>
+settings.
+
+<p>
+
+
+<dt><div class="code"><pre>
+g.region -up rast=soils
+</pre></div>
+
+<dd> The <b>-u</b> option suppresses the re-setting of the
+current region definition. This can be useful when it is
+desired to only extract region information. In this case,
+the cell header file for the soils map layer is printed
+without changing the current region settings.
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -up zoom=soils save=soils
+</pre></div>
+
+<dd> This will zoom into the smallest region which
+encompasses all non-NULL soils data values, and save the
+new region settings in a file to be called <em>soils</em>
+and stored under the <kbd>windows</kbd> directory in the
+user's current mapset. The current region settings are not
+changed.
+<p>
+
+
+<dt><div class="code"><pre>
+g.region b=0 t=3000 tbres=200 res3=100
+g.region -p3
+</pre></div>
+
+<dd> This will define the 3D region for voxel computations.
+In this example a volume with bottom (0m) to top (3000m)
+at horizontal resolution (100m) and vertical resolution (200m)
+is defined.
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -p
+</pre></div>
+
+<dd> This will print the current region in the format:
+
+<div class="code"><pre>
+projection: 1 (UTM)
+zone: 13
+datum: nad27
+ellipsoid: clark66
+north: 4928000
+south: 4914000
+west: 590000
+east: 609000
+nsres: 20
+ewres: 20
+rows: 700
+cols: 950
+</pre></div>
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -p3
+</pre></div>
+
+<dd> This will print the current region and the 3D region (used for voxels)
+in the format:
+
+<div class="code"><pre>
+projection: 1 (UTM)
+zone: 13
+datum: nad27
+ellipsoid: clark66
+north: 4928000
+south: 4914000
+west: 590000
+east: 609000
+top: 1.00000000
+bottom: 0.00000000
+nsres: 20
+nsres3: 20
+ewres: 20
+ewres3: 20
+tbres: 1
+rows: 700
+rows3: 700
+cols: 950
+cols3: 950
+depths: 1
+</pre></div>
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -g
+</pre></div>
+
+<dd> The <b>-g</b> option prints the region in the
+following script style (key=value) format:
+
+<div class="code"><pre>
+n=4928000
+s=4914000
+w=590000
+e=609000
+nsres=20
+ewres=20
+rows=700
+cols=950
+</pre></div>
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -bg
+</pre></div>
+
+<dd> The <b>-bg</b> option prints the region in the
+following script style (key=value) format plus the
+boundary box in latitude-longitude/WGS84:
+
+<div class="code"><pre>
+n=4928000
+s=4914000
+w=590000
+e=609000
+nsres=20
+ewres=20
+rows=700
+cols=950
+LL_W=-103.87080682
+LL_E=-103.62942884
+LL_N=44.50164277
+LL_S=44.37302019
+</pre></div>
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -l
+</pre></div>
+
+<dd> The <b>-l</b> option prints the region in the
+following format:
+
+<div class="code"><pre>
+long: -103.86789484 lat: 44.50165890 (north/west corner)
+long: -103.62895703 lat: 44.49904013 (north/east corner)
+long: -103.63190061 lat: 44.37303558 (south/east corner)
+long: -103.87032572 lat: 44.37564292 (south/west corner)
+rows: 700
+cols: 950
+Center longitude: 103:44:59.170374W [-103.74977]
+Center latitude: 44:26:14.439781N [44.43734]
+</pre></div>
+
+<p>
+
+<dt><div class="code"><pre>
+g.region -pm
+</pre></div>
+
+<dd> This will print the current region in the format
+ (latitude-longitude location):
+
+<div class="code"><pre>
+projection: 3 (Latitude-Longitude)
+zone: 0
+ellipsoid: wgs84
+north: 90N
+south: 40N
+west: 20W
+east: 20E
+nsres: 928.73944902
+ewres: 352.74269109
+rows: 6000
+cols: 4800
+</pre></div>
+Note that the resolution is here reported in meters, not decimal degrees.
+
+</dl>
+
+<p>
+Usage example of <em>g.region</em> in a shell with external software:<br>
+<!-- why not 'v.in.ogr spatial=' ?? -->
+Extract spatial subset of external vector map 'soils.shp' to new external
+vector map 'soils_cut.shp' using the OGR 'ogr2ogr' tool:<br>
+
+<div class="code"><pre>
+eval `g.region -g`
+ogr2ogr -spat $w $s $e $n soils_cut.shp soils.shp
+</pre></div>
+
+This requires that the location/SHAPE file projection match.
+
+
+<p>
+Usage example of <em>g.proj</em> and <em>g.region</em> in a shell with external software:<br>
+Extract spatial subset of external raster map 'p016r035_7t20020524_z17_nn30.tif'
+to new external raster map'p016r035_7t20020524_nc_spm_wake_nn30.tif using the GDAL
+'gdalwarp' tool:<br>
+
+<div class="code"><pre>
+eval `g.region -g`
+gdalwarp -t_srs "`g.proj -wf`" -te $w $s $e $n \
+ p016r035_7t20020524_z17_nn30.tif \
+ p016r035_7t20020524_nc_spm_wake_nn30.tif
+</pre></div>
+
+Here the input raster map does not have to match the location
+projection since it is reprojected on the fly.
+
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a HREF="d.zoom.html">d.zoom</a><br>
+<a HREF="g.access.html">g.access</a><br>
+<a HREF="g.mapsets.html">g.mapsets</a><br>
+<a HREF="g.proj.html">g.proj</a><BR>
+<a HREF="g.setproj.html">g.setproj</a><BR>
+environment variables: <a HREF="variables.html#internal">GRASS_REGION and WIND_OVERRIDE</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Michael Shapiro,
+U.S.Army Construction Engineering
+Research Laboratory
+
+<p>
+<i>Last changed: $Date$</i>
Copied: grass/trunk/general/g.region/local_proto.h (from rev 30698, grass/trunk/general/g.region/cmd/local_proto.h)
===================================================================
--- grass/trunk/general/g.region/local_proto.h (rev 0)
+++ grass/trunk/general/g.region/local_proto.h 2008-03-22 20:39:13 UTC (rev 30699)
@@ -0,0 +1,22 @@
+#ifndef GREGION_LOCAL_PROTO_H
+#define GREGION_LOCAL_PROTO_H
+
+#define PRINT_REG 0x01
+#define PRINT_SH 0x02
+#define PRINT_LL 0x04
+#define PRINT_EXTENT 0x08
+#define PRINT_CENTER 0x10
+#define PRINT_METERS 0x20
+#define PRINT_3D 0x40
+#define PRINT_MBBOX 0x80
+
+/* adjust.c */
+int adjust_window(struct Cell_head *, int, int, int);
+
+/* zoom.c */
+int zoom(struct Cell_head *, char *, char *);
+
+/* printwindow.c */
+int print_window(struct Cell_head *, int);
+
+#endif
Copied: grass/trunk/general/g.region/main.c (from rev 30698, grass/trunk/general/g.region/cmd/main.c)
===================================================================
--- grass/trunk/general/g.region/main.c (rev 0)
+++ grass/trunk/general/g.region/main.c 2008-03-22 20:39:13 UTC (rev 30699)
@@ -0,0 +1,891 @@
+/***************************************************************************
+ *
+ * MODULE: g.region (commandline)
+ * AUTHOR(S): Michael Shapiro, CERL
+ * datum added by Andreas Lange <andreas.lange at rhein-main.de>
+ * PURPOSE: Program to manage and print the boundary definitions for the
+ * geographic region.
+ *
+ * COPYRIGHT: (C) 2000 by the GRASS Development Team
+ *
+ * This program is free software under the GPL (>=v2)
+ * Read the file COPYING that comes with GRASS for details.
+ ****************************************************************************/
+
+#include <string.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/G3d.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+#include "projects.h"
+
+static int nsew(char *,char *,char *,char *);
+static void die(struct Option *);
+static char *llinfo(char *,char *,int);
+
+int main (int argc, char *argv[])
+{
+ int i;
+ int print_flag = 0;
+ int set_flag;
+ double x;
+ int ival;
+ int row_flag=0, col_flag=0;
+ struct Cell_head window, temp_window;
+ char *value;
+ char *name;
+ char *mapset;
+ char *err;
+ int projection;
+ char **rast_ptr, **vect_ptr;
+
+ struct GModule *module;
+ struct
+ {
+ struct Flag
+ *update,
+ *print,
+ *gprint,
+ *lprint,
+ *eprint,
+ *center,
+ *res_set,
+ *dist_res,
+ *dflt,
+ *z,
+ *savedefault,
+ *bbox;
+ } flag;
+ struct
+ {
+ struct Option
+ *north,*south,*east,*west,*top,*bottom,
+ *res, *nsres, *ewres, *res3, *tbres, *rows, *cols,
+ *save, *region, *view,
+ *raster, *raster3d, *align, *zoom, *vect;
+ } parm;
+
+ G_gisinit (argv[0]);
+
+ /* get current region.
+ * if current region not valid, set it from default
+ * note: G_get_default_window() dies upon error
+ */
+ G_get_default_window(&window);
+
+ module = G_define_module();
+ module->keywords = _("general");
+ module->description =
+ _("Manages the boundary definitions for the "
+ "geographic region.");
+
+ /* flags */
+
+ flag.dflt = G_define_flag();
+ flag.dflt->key = 'd';
+ flag.dflt->description = _("Set from default region");
+ flag.dflt->guisection = _("Existing");
+
+ flag.savedefault = G_define_flag();
+ flag.savedefault->key = 's';
+ flag.savedefault->description = _("Save as default region");
+ flag.savedefault->guisection = _("Existing");
+
+ flag.print = G_define_flag();
+ flag.print->key = 'p';
+ flag.print->description = _("Print the current region");
+ flag.print->guisection = _("Print");
+
+ flag.lprint = G_define_flag();
+ flag.lprint->key = 'l';
+ flag.lprint->description = _("Print the current region in lat/long "
+ "using the current ellipsoid/datum");
+ flag.lprint->guisection = _("Print");
+
+ flag.eprint = G_define_flag();
+ flag.eprint->key = 'e';
+ flag.eprint->description = _("Print the current region extent");
+ flag.eprint->guisection = _("Print");
+
+ flag.center = G_define_flag();
+ flag.center->key = 'c';
+ flag.center->description = _("Print the current region map center coordinates");
+ flag.center->guisection = _("Print");
+
+ flag.dist_res= G_define_flag();
+ flag.dist_res->key = 'm';
+ flag.dist_res->description = _("Print region resolution in meters (geodesic)");
+ flag.dist_res->guisection = _("Print");
+
+ flag.z = G_define_flag();
+ flag.z->key = '3';
+ flag.z->description = _("Print also 3D settings");
+ flag.z->guisection = _("Print");
+
+ flag.bbox = G_define_flag();
+ flag.bbox->key = 'b';
+ flag.bbox->description = _("Print the maximum bounding box in lat/long on WGS84");
+ flag.bbox->guisection = _("Print");
+
+ flag.gprint = G_define_flag();
+ flag.gprint->key = 'g';
+ flag.gprint->description = _("Print in shell script style");
+ flag.gprint->guisection = _("Print");
+
+ flag.res_set= G_define_flag();
+ flag.res_set->key = 'a';
+ flag.res_set->description = _("Align region to resolution (default = align to bounds, "
+ "works only for 2D resolution)");
+ flag.res_set->guisection = _("Bounds");
+
+ flag.update = G_define_flag();
+ flag.update->key = 'u';
+ flag.update->description = _("Do not update the current region");
+ flag.update->guisection = _("Effects");
+
+ /* parameters */
+
+ parm.region = G_define_option();
+ parm.region->key = "region";
+ parm.region->key_desc = "name";
+ parm.region->required = NO;
+ parm.region->multiple = NO;
+ parm.region->type = TYPE_STRING;
+ parm.region->description = _("Set current region from named region");
+ parm.region->gisprompt = "old,windows,region";
+ parm.region->guisection = _("Existing");
+
+ parm.raster = G_define_standard_option(G_OPT_R_MAP);
+ parm.raster->key = "rast";
+ parm.raster->required = NO;
+ parm.raster->multiple = YES;
+ parm.raster->description = _("Set region to match this raster map");
+ parm.raster->guisection = _("Existing");
+
+ parm.raster3d = G_define_standard_option(G_OPT_R3_MAP);
+ parm.raster3d->key = "rast3d";
+ parm.raster3d->required = NO;
+ parm.raster3d->multiple = NO;
+ parm.raster3d->description = _("Set region to match this 3D raster map (both 2D and 3D "
+ "values)");
+ parm.raster3d->guisection = _("Existing");
+
+ parm.vect = G_define_standard_option(G_OPT_V_MAP);
+ parm.vect->key = "vect";
+ parm.vect->required = NO;
+ parm.vect->multiple = YES;
+ parm.vect->description = _("Set region to match this vector map");
+ parm.vect->guisection = _("Existing");
+
+ parm.view = G_define_option();
+ parm.view->key = "3dview";
+ parm.view->key_desc = "name";
+ parm.view->required = NO;
+ parm.view->multiple = NO;
+ parm.view->type = TYPE_STRING;
+ parm.view->description = _("Set region to match this 3dview file");
+ parm.view->gisprompt = "old,3d.view,3d view";
+ parm.view->guisection = _("Existing");
+
+ parm.north = G_define_option();
+ parm.north->key = "n";
+ parm.north->key_desc = "value";
+ parm.north->required = NO;
+ parm.north->multiple = NO;
+ parm.north->type = TYPE_STRING;
+ parm.north->description = llinfo(_("Value for the northern edge"),
+ G_lat_format_string(), window.proj);
+ parm.north->guisection = _("Bounds");
+
+ parm.south = G_define_option();
+ parm.south->key = "s";
+ parm.south->key_desc = "value";
+ parm.south->required = NO;
+ parm.south->multiple = NO;
+ parm.south->type = TYPE_STRING;
+ parm.south->description = llinfo(_("Value for the southern edge"),
+ G_lat_format_string(), window.proj);
+ parm.south->guisection = _("Bounds");
+
+ parm.east = G_define_option();
+ parm.east->key = "e";
+ parm.east->key_desc = "value";
+ parm.east->required = NO;
+ parm.east->multiple = NO;
+ parm.east->type = TYPE_STRING;
+ parm.east->description = llinfo(_("Value for the eastern edge"),
+ G_lon_format_string(), window.proj);
+ parm.east->guisection = _("Bounds");
+
+ parm.west = G_define_option();
+ parm.west->key = "w";
+ parm.west->key_desc = "value";
+ parm.west->required = NO;
+ parm.west->multiple = NO;
+ parm.west->type = TYPE_STRING;
+ parm.west->description = llinfo(_("Value for the western edge"),
+ G_lon_format_string(), window.proj);
+ parm.west->guisection = _("Bounds");
+
+ parm.top = G_define_option();
+ parm.top->key = "t";
+ parm.top->key_desc = "value";
+ parm.top->required = NO;
+ parm.top->multiple = NO;
+ parm.top->type = TYPE_STRING;
+ parm.top->description = _("Value for the top edge");
+ parm.top->guisection = _("Bounds");
+
+ parm.bottom = G_define_option();
+ parm.bottom->key = "b";
+ parm.bottom->key_desc = "value";
+ parm.bottom->required = NO;
+ parm.bottom->multiple = NO;
+ parm.bottom->type = TYPE_STRING;
+ parm.bottom->description = _("Value for the bottom edge");
+ parm.bottom->guisection = _("Bounds");
+
+ parm.rows = G_define_option();
+ parm.rows->key = "rows";
+ parm.rows->key_desc = "value";
+ parm.rows->required = NO;
+ parm.rows->multiple = NO;
+ parm.rows->type = TYPE_INTEGER;
+ parm.rows->description = _("Number of rows in the new region");
+ parm.rows->guisection = _("Resolution");
+
+ parm.cols = G_define_option();
+ parm.cols->key = "cols";
+ parm.cols->key_desc = "value";
+ parm.cols->required = NO;
+ parm.cols->multiple = NO;
+ parm.cols->type = TYPE_INTEGER;
+ parm.cols->description = _("Number of columns in the new region");
+ parm.cols->guisection = _("Resolution");
+
+ parm.res = G_define_option();
+ parm.res->key = "res";
+ parm.res->key_desc = "value";
+ parm.res->required = NO;
+ parm.res->multiple = NO;
+ parm.res->type = TYPE_STRING;
+ parm.res->description = _("Grid resolution 2D (both north-south and east-west)");
+ parm.res->guisection = _("Resolution");
+
+ parm.res3 = G_define_option();
+ parm.res3->key = "res3";
+ parm.res3->key_desc = "value";
+ parm.res3->required = NO;
+ parm.res3->multiple = NO;
+ parm.res3->type = TYPE_STRING;
+ parm.res3->description = _("3D grid resolution (north-south, east-west and top-bottom)");
+ parm.res3->guisection = _("Resolution");
+
+ parm.nsres = G_define_option();
+ parm.nsres->key = "nsres";
+ parm.nsres->key_desc = "value";
+ parm.nsres->required = NO;
+ parm.nsres->multiple = NO;
+ parm.nsres->type = TYPE_STRING;
+ parm.nsres->description = llinfo(_("North-south grid resolution 2D"),
+ G_llres_format_string(), window.proj);
+ parm.nsres->guisection = _("Resolution");
+
+ parm.ewres = G_define_option();
+ parm.ewres->key = "ewres";
+ parm.ewres->key_desc = "value";
+ parm.ewres->required = NO;
+ parm.ewres->multiple = NO;
+ parm.ewres->type = TYPE_STRING;
+ parm.ewres->description = llinfo(_("East-west grid resolution 2D"),
+ G_llres_format_string(), window.proj);
+ parm.ewres->guisection = _("Resolution");
+
+ parm.tbres = G_define_option();
+ parm.tbres->key = "tbres";
+ parm.tbres->key_desc = "value";
+ parm.tbres->required = NO;
+ parm.tbres->multiple = NO;
+ parm.tbres->type = TYPE_STRING;
+ parm.tbres->description = _("Top-bottom grid resolution 3D");
+ parm.tbres->guisection = _("Resolution");
+
+ parm.zoom = G_define_option();
+ parm.zoom->key = "zoom";
+ parm.zoom->key_desc = "name";
+ parm.zoom->required = NO;
+ parm.zoom->multiple = NO;
+ parm.zoom->type = TYPE_STRING;
+ parm.zoom->description =
+ _("Shrink region until it meets non-NULL data from this raster map");
+ parm.zoom->gisprompt = "old,cell,raster";
+ parm.zoom->guisection = _("Bounds");
+
+ parm.align = G_define_option();
+ parm.align->key = "align";
+ parm.align->key_desc = "name";
+ parm.align->required = NO;
+ parm.align->multiple = NO;
+ parm.align->type = TYPE_STRING;
+ parm.align->description =
+ _("Adjust region cells to cleanly align with this raster map");
+ parm.align->gisprompt = "old,cell,raster";
+ parm.align->guisection = _("Bounds");
+
+ parm.save = G_define_option();
+ parm.save->key = "save";
+ parm.save->key_desc = "name";
+ parm.save->required = NO;
+ parm.save->multiple = NO;
+ parm.save->type = TYPE_STRING;
+ parm.save->description =
+ _("Save current region settings in named region file");
+ parm.save->gisprompt = "new,windows,region";
+ parm.save->guisection = _("Effects");
+
+ if (G_parser(argc,argv))
+ exit(EXIT_FAILURE);
+
+ projection = window.proj;
+
+ set_flag = ! flag.update->answer;
+
+ if (flag.print->answer)
+ print_flag |= PRINT_REG;
+
+ if (flag.gprint->answer)
+ print_flag |= PRINT_SH;
+
+ if (flag.lprint->answer)
+ print_flag |= PRINT_LL;
+
+ if (flag.eprint->answer)
+ print_flag |= PRINT_EXTENT;
+
+ if (flag.center->answer)
+ print_flag |= PRINT_CENTER;
+
+ if (flag.dist_res->answer)
+ print_flag |= PRINT_METERS;
+
+ if (flag.z->answer)
+ print_flag |= PRINT_3D;
+
+ if (flag.bbox->answer)
+ print_flag |= PRINT_MBBOX;
+
+ if (print_flag == PRINT_METERS)
+ print_flag |= PRINT_SH;
+
+ if (print_flag == PRINT_SH ||
+ print_flag & PRINT_3D ||
+ print_flag == PRINT_METERS + PRINT_SH)
+ {
+ print_flag |= PRINT_REG;
+ }
+
+ if (!flag.dflt->answer)
+ G_get_window(&window);
+
+ /* region= */
+ if ((name = parm.region->answer))
+ {
+ mapset = G_find_file ("windows", name, "");
+ if (!mapset)
+ G_fatal_error (_("Region <%s> not found"), name);
+ if (G__get_window (&window, "windows", name, mapset) != NULL)
+ G_fatal_error (_("Unable to read region <%s> in <%s>"), name, mapset);
+ }
+
+ /* 3dview= */
+ if ((name = parm.view->answer))
+ {
+ struct G_3dview v;
+ FILE *fp;
+ int ret;
+
+ mapset = G_find_file2 ("3d.view", name, "");
+ if (!mapset)
+ G_fatal_error (_("3dview file <%s> not found"), name);
+
+ G_3dview_warning(0); /* suppress boundary mismatch warning */
+
+ if(NULL == (fp = G_fopen_old("3d.view",name,mapset)))
+ G_fatal_error (_("Unable to open 3dview file <%s> in <%s>"), name, mapset);
+
+ G_copy (&temp_window, &window, sizeof(window));
+
+ if(0 > (ret = G_get_3dview(name, mapset, &v)))
+ G_fatal_error (_("Unable to read 3dview file <%s> in <%s>"), name, mapset);
+ if (ret == 0)
+ G_fatal_error (_("Old 3dview file. Region <%s> not found in <%s>"), name, mapset);
+
+
+ window.north = v.vwin.north;
+ window.south = v.vwin.south;
+ window.west = v.vwin.west;
+ window.east = v.vwin.east;
+
+ window.rows = v.vwin.rows;
+ window.cols = v.vwin.cols;
+ window.ns_res = v.vwin.ns_res;
+ window.ew_res = v.vwin.ew_res;
+
+ fclose (fp);
+
+ }
+
+ /* raster= */
+ if (parm.raster->answer)
+ {
+ int first = 0;
+ rast_ptr = parm.raster->answers;
+ for (; *rast_ptr != NULL; rast_ptr++)
+ {
+ char rast_name[GNAME_MAX];
+ strcpy (rast_name, *rast_ptr);
+ mapset = G_find_cell2 (rast_name, "");
+ if (!mapset)
+ G_fatal_error (_("Raster map <%s> not found"), rast_name);
+ if (G_get_cellhd (rast_name, mapset, &temp_window) < 0)
+ G_fatal_error (_("Unable to read header of raster map <%s@%s>"),
+ rast_name, mapset);
+ if (!first) {
+ G_copy (&window, &temp_window, sizeof(window));
+ first = 1;
+ } else {
+ window.north = (window.north > temp_window.north) ?
+ window.north : temp_window.north;
+ window.south = (window.south < temp_window.south) ?
+ window.south : temp_window.south;
+ window.east = (window.east > temp_window.east) ?
+ window.east : temp_window.east;
+ window.west = (window.west < temp_window.west) ?
+ window.west : temp_window.west;
+ }
+ }
+ G_adjust_Cell_head3(&window,0,0,0);
+ }
+
+
+ /* raster3d= */
+ if ((name = parm.raster3d->answer))
+ {
+ G3D_Region win;
+
+ if( (mapset = G_find_grid3(name, "")) == NULL )
+ G_fatal_error (_("3D raster map <%s> not found"), name);
+
+ if ( G3d_readRegionMap (name, mapset, &win) < 0 )
+ G_fatal_error (_("Unable to read header of 3D raster map <%s@%s>"),
+ name, mapset);
+
+ window.proj = win.proj;
+ window.zone = win.zone;
+ window.north = win.north;
+ window.south = win.south;
+ window.east = win.east;
+ window.west = win.west;
+ window.top = win.top;
+ window.bottom = win.bottom;
+ window.rows = win.rows;
+ window.rows3 = win.rows;
+ window.cols = win.cols;
+ window.cols3 = win.cols;
+ window.depths = win.depths;
+ window.ns_res = win.ns_res;
+ window.ns_res3 = win.ns_res;
+ window.ew_res = win.ew_res;
+ window.ew_res3 = win.ew_res;
+ window.tb_res = win.tb_res;
+ }
+
+ /* vect= */
+ if (parm.vect->answer)
+ {
+ int first = 0;
+ vect_ptr = parm.vect->answers;
+ for (; *vect_ptr != NULL; vect_ptr++)
+ {
+ struct Map_info Map;
+ BOUND_BOX box;
+ char vect_name[GNAME_MAX];
+ struct Cell_head map_window;
+
+ strcpy (vect_name, *vect_ptr);
+ mapset = G_find_vector2 (vect_name, "");
+ if (!mapset)
+ G_fatal_error (_("Vector map <%s> not found"), vect_name);
+
+ G_copy (&temp_window, &window, sizeof(window));
+
+ Vect_set_open_level (2);
+ if (2 != Vect_open_old (&Map, vect_name, mapset))
+ G_fatal_error (_("Unable to open vector map <%s@%s>"),
+ vect_name, mapset);
+
+ Vect_get_map_box (&Map, &box );
+ map_window = window;
+ map_window.north = box.N;
+ map_window.south = box.S;
+ map_window.west = box.W;
+ map_window.east = box.E;
+
+ if (!first)
+ {
+ G_copy (&window, &map_window, sizeof(window));
+ first = 1;
+ }
+ else
+ {
+ window.north = (window.north > map_window.north) ?
+ window.north : map_window.north;
+ window.south = (window.south < map_window.south) ?
+ window.south : map_window.south;
+ window.east = (window.east > map_window.east) ?
+ window.east : map_window.east;
+ window.west = (window.west < map_window.west) ?
+ window.west : map_window.west;
+ }
+
+ if(window.north == window.south)
+ {
+ window.north = window.north + 0.5 * temp_window.ns_res;
+ window.south = window.south - 0.5 * temp_window.ns_res;
+ }
+ if(window.east==window.west)
+ {
+ window.west = window.west - 0.5 * temp_window.ew_res;
+ window.east = window.east + 0.5 * temp_window.ew_res;
+ }
+
+ if(flag.res_set->answer)
+ G_align_window (&window, &temp_window);
+
+ Vect_close (&Map);
+ }
+ }
+
+ /* n= */
+ if ((value = parm.north->answer))
+ {
+ if((i = nsew(value, "n+", "n-", "s+")))
+ {
+ if (!G_scan_resolution (value+2, &x, window.proj))
+ die(parm.north);
+ switch(i)
+ {
+ case 1:
+ window.north += x;
+ break;
+ case 2:
+ window.north -= x;
+ break;
+ case 3:
+ window.north = window.south + x;
+ break;
+ }
+ }
+ else if (G_scan_northing (value, &x, window.proj))
+ window.north = x;
+ else
+ die(parm.north);
+ }
+
+ /* s= */
+ if ((value = parm.south->answer))
+ {
+ if((i = nsew(value, "s+", "s-", "n-")))
+ {
+ if (!G_scan_resolution (value+2, &x, window.proj))
+ die(parm.south);
+ switch(i)
+ {
+ case 1:
+ window.south += x;
+ break;
+ case 2:
+ window.south -= x;
+ break;
+ case 3:
+ window.south = window.north - x;
+ break;
+ }
+ }
+ else if (G_scan_northing (value, &x, window.proj))
+ window.south = x;
+ else
+ die(parm.south);
+ }
+
+ /* e= */
+ if ((value = parm.east->answer))
+ {
+ if((i = nsew(value, "e+", "e-", "w+")))
+ {
+ if (!G_scan_resolution (value+2, &x, window.proj))
+ die(parm.east);
+ switch(i)
+ {
+ case 1:
+ window.east += x;
+ break;
+ case 2:
+ window.east -= x;
+ break;
+ case 3:
+ window.east = window.west + x;
+ break;
+ }
+ }
+ else if (G_scan_easting (value, &x, window.proj))
+ window.east = x;
+ else
+ die(parm.east);
+ }
+
+ /* w= */
+ if ((value = parm.west->answer))
+ {
+ if((i = nsew(value, "w+", "w-", "e-")))
+ {
+ if (!G_scan_resolution (value+2, &x, window.proj))
+ die(parm.west);
+ switch(i)
+ {
+ case 1:
+ window.west += x;
+ break;
+ case 2:
+ window.west -= x;
+ break;
+ case 3:
+ window.west = window.east - x;
+ break;
+ }
+ }
+ else if (G_scan_easting (value, &x, window.proj))
+ window.west = x;
+ else
+ die(parm.west);
+ }
+
+ /* t= */
+ if ((value = parm.top->answer))
+ {
+ if((i = nsew(value, "t+", "t-", "b+")))
+ {
+ if ( sscanf (value+2, "%lf", &x) != 1 )
+ die(parm.top);
+ switch(i)
+ {
+ case 1:
+ window.top += x;
+ break;
+ case 2:
+ window.top -= x;
+ break;
+ case 3:
+ window.top = window.bottom + x;
+ break;
+ }
+ }
+ else if ( sscanf (value, "%lf", &x) == 1 )
+ window.top = x;
+ else
+ die(parm.top);
+ }
+
+ /* b= */
+ if ((value = parm.bottom->answer))
+ {
+ if((i = nsew(value, "b+", "b-", "t-")))
+ {
+ if ( sscanf (value+2, "%lf", &x) != 1 )
+ die(parm.bottom);
+ switch(i)
+ {
+ case 1:
+ window.bottom += x;
+ break;
+ case 2:
+ window.bottom -= x;
+ break;
+ case 3:
+ window.bottom = window.top - x;
+ break;
+ }
+ }
+ else if ( sscanf (value, "%lf", &x) == 1 )
+ window.bottom = x;
+ else
+ die(parm.bottom);
+ }
+
+ /* res= */
+ if ((value = parm.res->answer))
+ {
+ if (!G_scan_resolution (value, &x, window.proj))
+ die(parm.res);
+ window.ns_res = x;
+ window.ew_res = x;
+
+ if (flag.res_set->answer) {
+ window.north = ceil(window.north/x) * x ;
+ window.south = floor(window.south/x) * x ;
+ window.east = ceil(window.east/x) * x ;
+ window.west = floor(window.west/x) * x ;
+ }
+ }
+
+ /* res3= */
+ if ((value = parm.res3->answer))
+ {
+ if (!G_scan_resolution (value, &x, window.proj))
+ die(parm.res);
+ window.ns_res3 = x;
+ window.ew_res3 = x;
+ window.tb_res = x;
+ }
+
+ /* nsres= */
+ if ((value = parm.nsres->answer))
+ {
+ if (!G_scan_resolution (value, &x, window.proj))
+ die(parm.nsres);
+ window.ns_res = x;
+
+ if (flag.res_set->answer) {
+ window.north = ceil(window.north/x) * x ;
+ window.south = floor(window.south/x) * x ;
+ }
+ }
+
+ /* ewres= */
+ if ((value = parm.ewres->answer))
+ {
+ if (!G_scan_resolution (value, &x, window.proj))
+ die(parm.ewres);
+ window.ew_res = x;
+
+ if (flag.res_set->answer) {
+ window.east = ceil(window.east/x) * x ;
+ window.west = floor(window.west/x) * x ;
+ }
+ }
+
+ /* tbres= */
+ if ((value = parm.tbres->answer))
+ {
+ if ( sscanf (value, "%lf", &x) != 1 )
+ die(parm.tbres);
+ window.tb_res = x;
+
+ if (flag.res_set->answer) {
+ window.top = ceil(window.top/x) * x ;
+ window.bottom = floor(window.bottom/x) * x ;
+ }
+ }
+
+ /* rows= */
+ if ((value = parm.rows->answer))
+ {
+ if (sscanf (value, "%i", &ival) != 1 )
+ die(parm.rows);
+ window.rows = ival;
+ row_flag = 1;
+ }
+
+ /* cols= */
+ if ((value = parm.cols->answer))
+ {
+ if (sscanf (value, "%i", &ival) != 1 )
+ die(parm.cols);
+ window.cols = ival;
+ col_flag = 1;
+ }
+
+ /* zoom= */
+ if ((name = parm.zoom->answer))
+ {
+ mapset = G_find_cell2 (name, "");
+ if (!mapset)
+ G_fatal_error (_("Raster map <%s> not found"), name);
+ zoom (&window, name, mapset);
+ }
+
+ /* align= */
+ if ((name = parm.align->answer))
+ {
+ mapset = G_find_cell2 (name, "");
+ if (!mapset)
+ G_fatal_error (_("Raster map <%s> not found"), name);
+ if (G_get_cellhd (name, mapset, &temp_window) < 0)
+ G_fatal_error (_("Unable to read header of raster map <%s@%s>"),
+ name, mapset);
+ if ((err = G_align_window (&window, &temp_window)))
+ G_fatal_error (_("Raster map <%s@%s>: %s"), name, mapset, err);
+ }
+
+ /* save= */
+ if ((name = parm.save->answer))
+ {
+ if (G_legal_filename (name) < 0)
+ G_fatal_error (_("<%s> is an illegal region name"), name);
+ G_copy (&temp_window, &window, sizeof(window));
+ adjust_window (&temp_window,0,0,0);
+ if (G__put_window (&temp_window, "windows", name) < 0)
+ G_fatal_error (_("Unable to set region <%s>"), name);
+ }
+
+ adjust_window (&window,row_flag,col_flag,0);
+ if (set_flag)
+ {
+ if (G_put_window (&window) < 0)
+ G_fatal_error (_("Unable to update current region"));
+ }
+
+ if (flag.savedefault->answer) {
+ if (strcmp(G_mapset(),"PERMANENT") == 0) {
+ G__put_window( &window, "", "DEFAULT_WIND" );
+ }
+ else {
+ G_fatal_error(_("Unable to change default region. "
+ "The current mapset is not <PERMANENT>."));
+ }
+ } /* / flag.savedefault->answer */
+
+
+ if (print_flag)
+ print_window (&window, print_flag);
+
+ exit(EXIT_SUCCESS);
+}
+
+static void die(struct Option *parm)
+{
+ /*
+ G_usage();
+ */
+ G_fatal_error (_("Invalid input <%s=%s>"), parm->key, parm->answer);
+}
+
+static int nsew(char *value,char *a,char *b,char *c)
+{
+ if (strncmp (value, a, strlen(a)) == 0 ) return 1;
+ if (strncmp (value, b, strlen(b)) == 0 ) return 2;
+ if (strncmp (value, c, strlen(c)) == 0 ) return 3;
+ return 0;
+}
+
+static char *llinfo(char *msg,char *llformat,int proj)
+{
+ char buf[256];
+ if (proj != PROJECTION_LL)
+ return msg;
+
+ sprintf (buf, "%s (%s %s)", msg, _("format"), llformat);
+ return G_store(buf);
+}
Copied: grass/trunk/general/g.region/printwindow.c (from rev 30698, grass/trunk/general/g.region/cmd/printwindow.c)
===================================================================
--- grass/trunk/general/g.region/printwindow.c (rev 0)
+++ grass/trunk/general/g.region/printwindow.c 2008-03-22 20:39:13 UTC (rev 30699)
@@ -0,0 +1,667 @@
+#include <string.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/gprojects.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int print_window(struct Cell_head *window, int print_flag)
+{
+ char *prj, *datum, *ellps;
+ int x, width = 11;
+
+ char north[30], south[30], east[30], west[30], nsres[30], ewres[30],
+ nsres3[30], ewres3[30], tbres[30];
+ char buf[50];
+
+ double ew_dist1, ew_dist2, ns_dist1, ns_dist2;
+ double longitude, latitude;
+
+ if (print_flag & PRINT_SH)
+ x = -1;
+ else
+ x = window->proj;
+
+ G_format_northing(window->north, north, x);
+ G_format_northing(window->south, south, x);
+ G_format_easting(window->east, east, x);
+ G_format_easting(window->west, west, x);
+ G_format_resolution(window->ew_res, ewres, x);
+ G_format_resolution(window->ew_res3, ewres3, x);
+ G_format_resolution(window->ns_res, nsres, x);
+ G_format_resolution(window->ns_res3, nsres3, x);
+ G_format_resolution(window->tb_res, tbres, x);
+ G_begin_distance_calculations();
+
+ /* EW Dist at North edge */
+ ew_dist1 =
+ G_distance(window->east, window->north, window->west, window->north);
+ /* EW Dist at South Edge */
+ ew_dist2 =
+ G_distance(window->east, window->south, window->west, window->south);
+ /* NS Dist at East edge */
+ ns_dist1 =
+ G_distance(window->east, window->north, window->east, window->south);
+ /* NS Dist at West edge */
+ ns_dist2 =
+ G_distance(window->west, window->north, window->west, window->south);
+
+ /* width */
+ if (print_flag & PRINT_REG)
+ width = 11;
+
+ if (print_flag & PRINT_CENTER || print_flag & PRINT_MBBOX)
+ width = 16;
+
+ if (print_flag & PRINT_LL)
+ width = 18;
+
+ if (print_flag & PRINT_EXTENT)
+ width = 19;
+
+ /* flag.dist_res */
+ if (print_flag & PRINT_METERS)
+ {
+ sprintf(ewres, "%.8f", ((ew_dist1 + ew_dist2) / 2) / window->cols);
+ G_trim_decimal(ewres);
+ sprintf(ewres3, "%.8f", ((ew_dist1 + ew_dist2) / 2) / window->cols3);
+ G_trim_decimal(ewres3);
+ sprintf(nsres, "%.8f", ((ns_dist1 + ns_dist2) / 2) / window->rows);
+ G_trim_decimal(nsres);
+ sprintf(nsres3, "%.8f", ((ns_dist1 + ns_dist2) / 2) / window->rows3);
+ G_trim_decimal(nsres3);
+ sprintf(tbres, "%.8f", (window->top - window->bottom) / window->depths);
+ G_trim_decimal(tbres);
+ }
+
+ /* flag.print & flag.gprint */
+ if (print_flag & PRINT_REG) {
+ prj = G_database_projection_name();
+ if (!prj)
+ prj = "** unknown **";
+ /*
+ please remove before GRASS 7 released
+ backward compatibility issue
+
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "projection=%d\n", window->proj);
+ fprintf(stdout, "zone=%d\n", window->zone);
+ }
+ else
+ {
+ fprintf(stdout, "%-*s: %d (%s)\n", width, "projection", window->proj, prj);
+ fprintf(stdout, "%-*s: %d\n", width, "zone", window->zone);
+ }
+ */
+
+ if (! (print_flag & PRINT_SH))
+ {
+ fprintf(stdout, "%-*s %d (%s)\n", width, "projection:", window->proj, prj);
+ fprintf(stdout, "%-*s %d\n", width, "zone:", window->zone);
+ }
+
+ /* don't print datum/ellipsoid in XY-Locations */
+ if (window->proj != 0) {
+ datum = G_database_datum_name();
+ if (!datum)
+ datum = "** unknown (default: WGS84) **";
+ ellps = G_database_ellipse_name();
+ if (!ellps)
+ ellps = "** unknown (default: WGS84) **";
+
+ /*
+ please remove before GRASS 7 released
+ backward compatibility issue
+
+ if (print_flag & PRINT_SH)
+ {
+ if (datum[0] != '*')
+ fprintf(stdout, "datum=%s\n", datum);
+ else
+ fprintf(stdout, "datum=wgs84\n");
+ if (ellps[0] != '*')
+ fprintf(stdout, "ellipsoid=%s\n", ellps);
+ else
+ fprintf(stdout, "ellipsoid=wgs84\n");
+ }
+ else
+ {
+ fprintf(stdout, "%-*s %s\n", width, "datum:", datum);
+ fprintf(stdout, "%-*s %s\n", width, "ellipsoid:", ellps);
+ }
+ */
+
+ if (! (print_flag & PRINT_SH))
+ {
+ fprintf(stdout, "%-*s %s\n", width, "datum:", datum);
+ fprintf(stdout, "%-*s %s\n", width, "ellipsoid:", ellps);
+ }
+ }
+
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "n=%s\n", north);
+ fprintf(stdout, "s=%s\n", south);
+ fprintf(stdout, "w=%s\n", west);
+ fprintf(stdout, "e=%s\n", east);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "t=%g\n", window->top);
+ fprintf(stdout, "b=%g\n", window->bottom);
+ }
+ fprintf(stdout, "nsres=%s\n", nsres);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "nsres3=%s\n", nsres3);
+ }
+ fprintf(stdout, "ewres=%s\n", ewres);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "ewres3=%s\n", ewres3);
+ fprintf(stdout, "tbres=%s\n", tbres);
+ }
+ fprintf(stdout, "rows=%d\n", window->rows);
+ if (print_flag & PRINT_3D)
+ fprintf(stdout, "rows3=%d\n", window->rows3);
+ fprintf(stdout, "cols=%d\n", window->cols);
+ if (print_flag & PRINT_3D)
+ {
+ fprintf(stdout, "cols3=%d\n", window->cols3);
+ fprintf(stdout, "depths=%d\n", window->depths);
+ }
+#ifdef HAVE_LONG_LONG_INT
+ fprintf(stdout, "cells=%lld\n", (long long) window->rows * window->cols);
+ if (print_flag & PRINT_3D)
+ fprintf(stdout, "3dcells=%lld\n",
+ (long long) window->rows3 * window->cols3 * window -> depths);
+#else
+ fprintf(stdout, "cells=%ld\n", (long) window->rows * window->cols);
+ if (print_flag & PRINT_3D)
+ fprintf(stdout, "3dcells=%ld\n",
+ (long) window->rows3 * window->cols3 * window -> depths);
+#endif
+ }
+ else
+ {
+ fprintf(stdout, "%-*s %s\n", width, "north:", north);
+ fprintf(stdout, "%-*s %s\n", width, "south:", south);
+ fprintf(stdout, "%-*s %s\n", width, "west:", west);
+ fprintf(stdout, "%-*s %s\n", width, "east:", east);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "%-*s %.8f\n", width, "top:", window->top);
+ fprintf(stdout, "%-*s %.8f\n", width, "bottom:", window->bottom);
+ }
+ fprintf(stdout, "%-*s %s\n", width, "nsres:", nsres);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "%-*s %s\n", width, "nsres3:", nsres3);
+ }
+ fprintf(stdout, "%-*s %s\n", width, "ewres:", ewres);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "%-*s %s\n", width, "ewres3:", ewres3);
+ fprintf(stdout, "%-*s %s\n", width, "tbres:", tbres);
+ }
+
+ fprintf(stdout, "%-*s %d\n", width, "rows:", window->rows);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "%-*s %d\n", width, "rows3:", window->rows3);
+ }
+ fprintf(stdout, "%-*s %d\n", width, "cols:", window->cols);
+ if (print_flag & PRINT_3D) {
+ fprintf(stdout, "%-*s %d\n", width, "cols3:", window->cols3);
+ fprintf(stdout, "%-*s %d\n", width, "depths:", window->depths);
+ }
+#ifdef HAVE_LONG_LONG_INT
+ fprintf(stdout, "%-*s %lld\n", width, "cells:",
+ (long long) window->rows * window->cols);
+ if (print_flag & PRINT_3D)
+ fprintf(stdout, "%-*s %lld\n", width, "3dcells:",
+ (long long) window->rows3 * window->cols3 * window -> depths);
+#else
+ fprintf(stdout, "%-*s %ld\n", width, "cells:",
+ (long) window->rows * window->cols);
+ if (print_flag & PRINT_3D)
+ fprintf(stdout, "%-*s %ld\n", width, "3dcells:",
+ (long) window->rows3 * window->cols3 * window -> depths);
+#endif
+ }
+ }
+
+ /* flag.lprint: show boundaries in lat/long MN 2001 */
+ if (print_flag & PRINT_LL)
+ {
+ double lo1, la1, lo2, la2, lo3, la3, lo4, la4;
+ double mid_n_lo, mid_n_la, mid_s_lo, mid_s_la, mid_w_lo, mid_w_la,
+ mid_e_lo, mid_e_la;
+
+ /* if coordinates are not in lat/long format, transform them: */
+ if ((G_projection() != PROJECTION_LL) && window->proj != 0) {
+ /* projection information of input map */
+ struct Key_Value *in_proj_info, *in_unit_info;
+ struct pj_info iproj; /* input map proj parameters */
+ struct pj_info oproj; /* output map proj parameters */
+
+ /* 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"));
+
+ G_free_key_value(in_proj_info);
+ G_free_key_value(in_unit_info);
+
+ /* 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.;
+ mid_w_lo = (lo1 + lo4) / 2.; /* not needed */
+ mid_w_la = (la1 + la4) / 2.; /* not needed */
+ mid_e_lo = (lo3 + lo2) / 2.; /* not needed */
+ mid_e_la = (la3 + la2) / 2.; /* not needed */
+
+ G_debug(3, "mid_n_lo %f", mid_n_lo);
+ G_debug(3, "mid_s_lo %f", mid_s_lo);
+ G_debug(3, "mid_n_la %f", mid_n_la);
+ G_debug(3, "mid_s_la %f", mid_s_la);
+ G_debug(3, "mid_w_lo %f", mid_w_lo);
+ G_debug(3, "mid_e_lo %f", mid_e_lo);
+ G_debug(3, "mid_w_la %f", mid_w_la);
+ G_debug(3, "mid_e_la %f", mid_e_la);
+
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "nw_long=%.8f\nnw_lat=%.8f\n", lo1, la1);
+ fprintf(stdout, "ne_long=%.8f\nne_lat=%.8f\n", lo2, la2);
+ fprintf(stdout, "se_long=%.8f\nse_lat=%.8f\n", lo3, la3);
+ fprintf(stdout, "sw_long=%.8f\nsw_lat=%.8f\n", lo4, la4);
+ G_format_easting((mid_n_lo + mid_s_lo) / 2., buf, PROJECTION_LL);
+ fprintf(stdout, "center_long=%.8f\n",
+ (mid_n_lo + mid_s_lo) / 2.);
+ G_format_northing((mid_n_la + mid_s_la) / 2., buf, PROJECTION_LL);
+ fprintf(stdout, "center_lat=%.8f\n",
+ (mid_n_la + mid_s_la) / 2.);
+
+ }
+ else
+ {
+ G_format_easting(lo1, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s long: %s ", width, "north-west corner:", buf);
+ G_format_northing(la1, buf, PROJECTION_LL);
+ fprintf(stdout, "lat: %s\n", buf);
+
+ G_format_easting(lo2, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s long: %s ", width, "north-east corner:", buf);
+ G_format_northing(la2, buf, PROJECTION_LL);
+ fprintf(stdout, "lat: %s\n", buf);
+
+ G_format_easting(lo3, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s long: %s ", width, "south-east corner:", buf);
+ G_format_northing(la3, buf, PROJECTION_LL);
+ fprintf(stdout, "lat: %s\n", buf);
+
+ G_format_easting(lo4, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s long: %s ", width, "south-west corner:", buf);
+ G_format_northing(la4, buf, PROJECTION_LL);
+ fprintf(stdout, "lat: %s\n", buf);
+
+ G_format_easting((mid_n_lo + mid_s_lo) / 2., buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %11s\n", width, "center longitude:", buf);
+ G_format_northing((mid_n_la + mid_s_la) / 2., buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %11s\n", width, "center latitude:", buf);
+ }
+
+ if (!(print_flag & PRINT_REG))
+ {
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "rows=%d\n", window->rows);
+ fprintf(stdout, "cols=%d\n", window->cols);
+ }
+ else
+ {
+ fprintf(stdout, "%-*s %d\n", width, "rows:", window->rows);
+ fprintf(stdout, "%-*s %d\n", width, "cols:", window->cols);
+ }
+ }
+ }
+ else /* in lat/long already */
+ {
+ if (window->proj != 0)
+ G_message (_("You are already in lat/long. Use -p flag instead."));
+ else
+ G_message (_("You are in xy location (no projection possible, "
+ "use -p flag instead)."));
+ }
+ }
+
+ /* flag.eprint */
+ if (print_flag & PRINT_EXTENT)
+ {
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "ns_extent=%f\n", window->north - window->south);
+ fprintf(stdout, "ew_extent=%f\n", window->east - window->west);
+ }
+ else
+ {
+ if (G_projection() != PROJECTION_LL)
+ {
+ fprintf(stdout, "%-*s %f\n", width, "north-south extent:",
+ window->north - window->south);
+ fprintf(stdout, "%-*s %f\n", width, "east-west extent:",
+ window->east - window->west);
+ }
+ else
+ {
+ G_format_northing(window->north - window->south, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "north-south extent:",
+ buf);
+ G_format_easting(window->east - window->west, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "east-west extent:",
+ buf);
+ }
+ }
+ }
+
+ /* flag.center */
+ if (print_flag & PRINT_CENTER)
+ {
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "center_easting=%f\n",
+ ((window->west - window->east) / 2. + window->east));
+ fprintf(stdout, "center_northing=%f\n",
+ ((window->north - window->south) / 2. + window->south));
+ }
+ else
+ {
+ if (G_projection() != PROJECTION_LL)
+ {
+ fprintf(stdout, "%-*s %f\n", width, "center easting:",
+ ((window->west - window->east) / 2. + window->east));
+ fprintf(stdout, "%-*s %f\n", width, "center northing:",
+ ((window->north - window->south) / 2. + window->south));
+ }
+ else
+ {
+ G_format_northing((window->north - window->south) / 2. + window->south,
+ buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "north-south center:", buf);
+ G_format_easting((window->west - window->east) / 2. + window->east,
+ buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "east-west center:", buf);
+ }
+ }
+ }
+
+ /* flag.bbox
+ Calculate the largest boudingbox in lat-lon coordinates
+ and print it to stdout
+ */
+ if (print_flag & PRINT_MBBOX)
+ {
+ double lo1, la1, lo2, la2, lo3, la3, lo4, la4;
+ double sh_ll_w, sh_ll_e, sh_ll_n, sh_ll_s;
+ /*double sh_ll_rows, sh_ll_cols; */
+
+ /* Needed to calculate the LL bounding box */
+ if ((G_projection() != PROJECTION_XY)) {
+ /* projection information of input and output map */
+ struct Key_Value *in_proj_info, *in_unit_info, *out_proj_info, *out_unit_info;
+ struct pj_info iproj; /* input map proj parameters */
+ struct pj_info oproj; /* output map proj parameters */
+ char buff[100], dum[100];
+
+ /* 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 and wgs84 ellipsoid */
+ out_proj_info = G_create_key_value();
+ out_unit_info = G_create_key_value();
+
+ G_set_key_value("proj", "ll", out_proj_info);
+
+ if (G_get_datumparams_from_projinfo(in_proj_info, buff, dum) < 0)
+ G_fatal_error(_
+ ("WGS84 output not possible as this location does not contain "
+ "datum transformation parameters. Try running g.setproj."));
+ else
+ G_set_key_value("datum", "wgs84", out_proj_info);
+
+ G_set_key_value("unit", "degree", out_unit_info);
+ G_set_key_value("units", "degrees", out_unit_info);
+ G_set_key_value("meters", "1.0", out_unit_info);
+
+ if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
+ G_fatal_error(_
+ ("Unable to update lat/long projection parameters"));
+
+ G_free_key_value(in_proj_info);
+ G_free_key_value(in_unit_info);
+ G_free_key_value(out_proj_info);
+ G_free_key_value(out_unit_info);
+
+
+ /* 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 up lat/long projection parameters"));
+ */
+
+ /* 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;
+
+ /*Calculate the largest bounding box */
+ if (fabs(lo3) > fabs(lo4)) {
+ if (fabs(lo4) < fabs(lo1))
+ sh_ll_w = lo4;
+ else
+ sh_ll_w = lo1;
+ if (fabs(lo3) > fabs(lo2))
+ sh_ll_e = lo3;
+ else
+ sh_ll_e = lo2;
+ }
+ else {
+ if (fabs(lo4) > fabs(lo1))
+ sh_ll_w = lo4;
+ else
+ sh_ll_w = lo1;
+ if (fabs(lo3) < fabs(lo2))
+ sh_ll_e = lo3;
+ else
+ sh_ll_e = lo2;
+ }
+
+ if (fabs(la4) < fabs(la1)) {
+ if (fabs(la1) > fabs(la2))
+ sh_ll_n = la1;
+ else
+ sh_ll_n = la2;
+ if (fabs(la4) < fabs(la3))
+ sh_ll_s = la4;
+ else
+ sh_ll_s = la3;
+ }
+ else {
+ if (fabs(la1) < fabs(la2))
+ sh_ll_n = la1;
+ else
+ sh_ll_n = la2;
+ if (fabs(la4) > fabs(la3))
+ sh_ll_s = la4;
+ else
+ sh_ll_s = la3;
+ }
+
+ /* print the larg bounding box */
+ if (print_flag & PRINT_SH)
+ {
+ fprintf(stdout, "ll_n=%.8f\n", sh_ll_n);
+ fprintf(stdout, "ll_s=%.8f\n", sh_ll_s);
+ fprintf(stdout, "ll_w=%.8f\n", sh_ll_w);
+ fprintf(stdout, "ll_e=%.8f\n", sh_ll_e);
+ fprintf(stdout, "ll_clon=%.8f\n", (sh_ll_e - sh_ll_w)/2. + sh_ll_w);
+ fprintf(stdout, "ll_clat=%.8f\n", (sh_ll_n - sh_ll_s)/2. + sh_ll_s);
+ }
+ else
+ {
+ G_format_northing(sh_ll_n, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "north longitude:", buf);
+ G_format_northing(sh_ll_s, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "south longitude:", buf);
+ G_format_easting(sh_ll_w, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "west latitude:", buf);
+ G_format_easting(sh_ll_e, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "east latitude:", buf);
+ G_format_easting((sh_ll_e - sh_ll_w)/2. + sh_ll_w, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "center latitude:", buf);
+ G_format_northing((sh_ll_n - sh_ll_s)/2. + sh_ll_s, buf, PROJECTION_LL);
+ fprintf(stdout, "%-*s %s\n", width, "center longitude:", buf);
+ }
+
+ /*It should be calculated which number of rows and cols we have in LL */
+ /*
+ fprintf (stdout, "LL_ROWS=%f \n",sh_ll_rows);
+ fprintf (stdout, "LL_COLS=%f \n",sh_ll_cols);
+ */
+ }
+ else
+ {
+ G_warning (_("No LatLong information for XY-Projection"));
+ }
+ }
+
+ return 0;
+}
Copied: grass/trunk/general/g.region/zoom.c (from rev 30698, grass/trunk/general/g.region/cmd/zoom.c)
===================================================================
--- grass/trunk/general/g.region/zoom.c (rev 0)
+++ grass/trunk/general/g.region/zoom.c 2008-03-22 20:39:13 UTC (rev 30699)
@@ -0,0 +1,75 @@
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "local_proto.h"
+
+int zoom (struct Cell_head *window,char *name,char *mapset)
+{
+ int fd;
+ void *raster, *rast_ptr;
+ RASTER_MAP_TYPE map_type;
+ int row, col;
+ int nrows, ncols;
+ int top,bottom,left,right,mark;
+ double north, south, east, west;
+
+ adjust_window (window,0,0,0);
+ G_set_window (window);
+ nrows = window->rows;
+ ncols = window->cols;
+
+ fd = G_open_cell_old (name, mapset);
+ if (fd < 0)
+ G_fatal_error (_("Unable to open raster map <%s> in <%s>"),
+ name, mapset);
+ map_type = G_get_raster_map_type(fd);
+ raster = G_allocate_raster_buf(map_type);
+
+/* find first non-null row */
+ top = nrows;
+ bottom = -1;
+ left = ncols;
+ right = -1;
+ for (row = 0; row < nrows; row++)
+ {
+ if (G_get_raster_row (fd, rast_ptr = raster,row,map_type) < 0)
+ G_fatal_error (_("Could not read from <%s>"),name);
+ for (col = 0; col < ncols; col++)
+ {
+ if (!G_is_null_value(rast_ptr, map_type))
+ break;
+ rast_ptr = G_incr_void_ptr(rast_ptr, G_raster_size(map_type));
+ }
+ if (col == ncols)
+ continue;
+ if (row < top) top = row;
+ if (row > bottom) bottom = row;
+ if (col < left) left = col;
+ for (mark = col; col < ncols; col++)
+ {
+ if (!G_is_null_value(rast_ptr, map_type))
+ mark = col;
+ rast_ptr = G_incr_void_ptr(rast_ptr, G_raster_size(map_type));
+ }
+ if (mark > right)
+ right = mark;
+ }
+ G_close_cell (fd);
+ G_free (raster);
+
+/* no data everywhere? */
+ if (bottom < 0)
+ return 0;
+
+ north = window->north - top * window->ns_res;
+ south = window->north - (bottom+1) * window->ns_res;
+ west = window->west + left * window->ew_res;
+ east = window->west + (right+1) * window->ew_res;
+
+ window->north = north;
+ window->south = south;
+ window->east = east ;
+ window->west = west ;
+
+ return 1;
+}
More information about the grass-commit
mailing list