[GRASS-SVN] r70778 - in grass-addons/grass6/raster: . r.fill.gaps
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Mar 21 04:26:32 PDT 2017
Author: benducke
Date: 2017-03-21 04:26:32 -0700 (Tue, 21 Mar 2017)
New Revision: 70778
Added:
grass-addons/grass6/raster/r.fill.gaps/
grass-addons/grass6/raster/r.fill.gaps/Makefile
grass-addons/grass6/raster/r.fill.gaps/cell_funcs.c
grass-addons/grass6/raster/r.fill.gaps/cell_funcs.h
grass-addons/grass6/raster/r.fill.gaps/description.html
grass-addons/grass6/raster/r.fill.gaps/main.c
grass-addons/grass6/raster/r.fill.gaps/r.fill.gaps.01.png
grass-addons/grass6/raster/r.fill.gaps/r.fill.gaps.02.png
Log:
Initial commit of r.fill.gaps:
Quickly fills small gaps in dense, rasterized data.
Added: grass-addons/grass6/raster/r.fill.gaps/Makefile
===================================================================
--- grass-addons/grass6/raster/r.fill.gaps/Makefile (rev 0)
+++ grass-addons/grass6/raster/r.fill.gaps/Makefile 2017-03-21 11:26:32 UTC (rev 70778)
@@ -0,0 +1,11 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.fill.gaps
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Added: grass-addons/grass6/raster/r.fill.gaps/cell_funcs.c
===================================================================
--- grass-addons/grass6/raster/r.fill.gaps/cell_funcs.c (rev 0)
+++ grass-addons/grass6/raster/r.fill.gaps/cell_funcs.c 2017-03-21 11:26:32 UTC (rev 70778)
@@ -0,0 +1,140 @@
+/***************************************************************************
+ *
+ * MODULE: r.fill.gaps (commandline)
+ * FILE: cell_funcs.c
+ * AUTHOR(S): Benjamin Ducke
+ * PURPOSE: Provide cell type dependent functions for cell data operations.
+ * At program initialization, a function pointer is set to point to
+ * the right one of the three alternatives.
+ * This function pointer can later be used to work with cell data of
+ * different types without having to go through any switching logics.
+ *
+ * COPYRIGHT: (C) 2011 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 <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "cell_funcs.h"
+
+
+/*
+ * Write cell values.
+ */
+void write_cell_value_c ( void *cell_output, void *cell_input ) {
+ G_set_raster_value_c(cell_output,G_get_raster_value_c(cell_input,IN_TYPE), OUT_TYPE);
+}
+
+void write_cell_value_f ( void *cell_output, void *cell_input ) {
+ G_set_raster_value_f(cell_output,G_get_raster_value_f(cell_input,IN_TYPE), OUT_TYPE);
+}
+
+void write_cell_value_d ( void *cell_output, void *cell_input ) {
+ G_set_raster_value_d(cell_output,G_get_raster_value_d(cell_input,IN_TYPE), OUT_TYPE);
+}
+
+
+/*
+ * Write a double value into a cell (truncates for CELL type output).
+ */
+
+void write_double_value_c ( void *cell, double val ) {
+ G_set_raster_value_c(cell,(CELL)val, OUT_TYPE);
+}
+
+void write_double_value_f ( void *cell, double val ) {
+ G_set_raster_value_f(cell,(FCELL)val, OUT_TYPE);
+}
+
+void write_double_value_d ( void *cell, double val ) {
+ G_set_raster_value_d(cell,(DCELL)val, OUT_TYPE);
+}
+
+
+
+/*
+ * Check for "no data".
+ */
+int is_null_value_c ( void *cell ) {
+ return (G_is_c_null_value ((CELL*)cell));
+}
+
+int is_null_value_f ( void *cell ) {
+ return (G_is_f_null_value ((FCELL*)cell));
+}
+
+int is_null_value_d ( void *cell ) {
+ return (G_is_d_null_value ((DCELL*)cell));
+}
+
+
+/*
+ * Set consecutive cells to "no data".
+ */
+void set_null_c (void* cells, unsigned long count) {
+ G_set_c_null_value ((CELL*)cells, count);
+}
+
+void set_null_f (void* cells, unsigned long count) {
+ G_set_f_null_value ((FCELL*)cells, count);
+}
+
+void set_null_d (void* cells, unsigned long count) {
+ G_set_d_null_value ((DCELL*)cells, count);
+}
+
+
+/*
+ * Call this init function once it is clear which input
+ * and output type will be used (CELL, DCELL or FCELL).
+ * I.e. right after IN_TYPE and OUT_TYPE have been set
+ * to valid values;
+ */
+void init_cell_funcs () {
+
+ CELL_ERR_SIZE = sizeof (FCELL);
+
+ if ( IN_TYPE == CELL_TYPE ) {
+ CELL_IN_SIZE = sizeof (CELL);
+ CELL_IN_PTR_SIZE = sizeof (CELL*);
+ IS_NULL = &is_null_value_c;
+ }
+ if ( IN_TYPE == FCELL_TYPE ) {
+ CELL_IN_SIZE = sizeof (FCELL);
+ CELL_IN_PTR_SIZE = sizeof (FCELL*);
+ IS_NULL = &is_null_value_f;
+ }
+ if ( IN_TYPE == DCELL_TYPE ) {
+ CELL_IN_SIZE = sizeof (DCELL);
+ CELL_IN_PTR_SIZE = sizeof (DCELL*);
+ IS_NULL = &is_null_value_d;
+ }
+ if ( OUT_TYPE == CELL_TYPE ) {
+ CELL_OUT_SIZE = sizeof (CELL);
+ CELL_OUT_PTR_SIZE = sizeof (CELL*);
+ WRITE_CELL_VAL = &write_cell_value_c;
+ WRITE_DOUBLE_VAL = &write_double_value_c;
+ SET_NULL = &set_null_c;
+ }
+ if ( OUT_TYPE == FCELL_TYPE ) {
+ CELL_OUT_SIZE = sizeof (FCELL);
+ CELL_OUT_PTR_SIZE = sizeof (FCELL*);
+ WRITE_CELL_VAL = &write_cell_value_f;
+ WRITE_DOUBLE_VAL = &write_double_value_f;
+ SET_NULL = &set_null_f;
+ }
+ if ( OUT_TYPE == DCELL_TYPE ) {
+ CELL_OUT_SIZE = sizeof (DCELL);
+ CELL_OUT_PTR_SIZE = sizeof (DCELL*);
+ WRITE_CELL_VAL = &write_cell_value_d;
+ WRITE_DOUBLE_VAL = &write_double_value_d;
+ SET_NULL = &set_null_d;
+ }
+}
Added: grass-addons/grass6/raster/r.fill.gaps/cell_funcs.h
===================================================================
--- grass-addons/grass6/raster/r.fill.gaps/cell_funcs.h (rev 0)
+++ grass-addons/grass6/raster/r.fill.gaps/cell_funcs.h 2017-03-21 11:26:32 UTC (rev 70778)
@@ -0,0 +1,32 @@
+/***************************************************************************
+ *
+ * MODULE: r.fill.gaps (commandline)
+ * FILE: cell_funcs.h
+ * AUTHOR(S): Benjamin Ducke
+ * PURPOSE: See cell_funcs.c
+ *
+ * COPYRIGHT: (C) 2011 by the GRASS Development Team
+ *
+ * This program is free software under the GPL (>=v2)
+ * Read the file COPYING that comes with GRASS for details.
+ ****************************************************************************
+ */
+
+#ifndef CELL_FUNCS_H
+#define CELL_FUNCS_H
+
+RASTER_MAP_TYPE IN_TYPE; /* stuff for cell input and output data */
+RASTER_MAP_TYPE OUT_TYPE;
+/* below are the sizes (in bytes) for the different GRASS cell types */
+unsigned char CELL_IN_SIZE;
+unsigned char CELL_IN_PTR_SIZE;
+unsigned char CELL_OUT_SIZE;
+unsigned char CELL_OUT_PTR_SIZE;
+unsigned char CELL_ERR_SIZE;
+
+void (*WRITE_CELL_VAL) (void*, void*); /* write a cell value of any type into an output cell */
+void (*WRITE_DOUBLE_VAL) (void*, double); /* write a double value into an output cell */
+int (*IS_NULL) (void*); /* check if a cell is "no data" */
+void (*SET_NULL) (void*, unsigned long); /* write null value(s) */
+
+#endif /* CELL_FUNCS_H */
Added: grass-addons/grass6/raster/r.fill.gaps/description.html
===================================================================
--- grass-addons/grass6/raster/r.fill.gaps/description.html (rev 0)
+++ grass-addons/grass6/raster/r.fill.gaps/description.html 2017-03-21 11:26:32 UTC (rev 70778)
@@ -0,0 +1,255 @@
+<h2>DESCRIPTION</h2>
+
+
+<em><b>r.fill.gaps</b></em> - Fast gap filling and interpolation of dense raster data.
+
+<p>
+
+The <em>r.fill.gaps</em> module is capable of quickly filling small "no data" areas (gaps) in large raster maps.
+As a rule of thumb, there should be at least as many data cells as there are "no data" cells in
+the input raster map. Several interpolation modes are available.
+By default, all values of the input raster map will be replaced with locally interpolated values in the output raster map.
+With dense data and small gaps, this should result in only slight deviations from the original data
+and produce smooth output. Alternatively, setting the <em>-p</em> flag will preserve the original cell values.
+
+<p>
+
+Available gap filling modes:
+<ul>
+<li>spatially weighted mean (default)</li>
+<li>simple mean</li>
+<li>simple median</li>
+<li>simple mode</li>
+</ul>
+
+<p>
+
+The spatially weighted mean is equivalent to an Inverse Distance Weighting (IDW; see
+also <em><a href="r.surf.idw.html">r.surf.idw</a></em>) interpolation. The simple
+mean is equivalent to a low-pass filter. Median and mode replacements can also be
+achieved using <em><a href="r.neighbors.html">r.neighbors</a></em>.
+
+<p>
+
+Note that <em>r.fill.gaps</em> is highly optimized for fast processing of large raster datasets
+with a small interpolation distance threshold (i.e. closing small gaps). As a trade-off
+for speed and a small memory foot print, some spatial accuracy is lost due to the rounding
+of the distance threshold to the closest approximation in input raster cells and the use of a matrix of precomputed weights at cell resolution (see further
+down for details). Note also that processing time will increase exponentially with higher distance settings.
+Cells outside the distance threshold will not be interpolated, preserving the edges of the data zones in the input data.
+
+<p>
+
+This module is not well suited for interpolating sparse input data and closing large gaps.
+For such purposes solutions more appropriate methods are available, such as <em><a href="v.surf.idw.html">v.surf.idw</a></em>
+or <em><a href="v.surf.rst.html">v.surf.rst</a></em>.
+
+<p>
+
+Applications where the properties of <em>r.fill.gaps</em> are advantageous include the processing of
+high resolution, close range scanning and remote sensing datasets. Such datasets commonly feature densely spaced
+measurements that have some gaps after rasterization, due to blind spots, instrument failures,
+and misalignments between the GIS' raster cell grids and the original measurement locations.
+In these cases, <em>r.fill.gaps</em> should typically be run using the "weighted mean" (default) method
+and with a small distance setting (the default is to use a search radius of just three cells).
+
+<p>
+
+The images below show a gradiometer dataset with gaps and its interpolated equivalent, produced using
+the spatially weighted mean operator (<em>mode="wmean"</em>).
+
+<p>
+
+<img src="r.fill.gaps.01.png" alt=""> <img src="r.fill.gaps.02.png" alt="">
+
+<p>
+
+In addition, <em>r.fill.gaps</em> can be useful for raster map generalization. Often, this involves removing small
+clumps of categorized cells and then filling the resulting data gaps without "inventing" any new values.
+In such cases, the "mode" or "median" interpolators should be used.
+
+
+<h2>USAGE</h2>
+
+The most critical user-provided settings are the interpolation/gap filling method (<em>mode=</em>)
+and the maximum distance, up to which <em>r.fill.gaps</em> will scan for given data points
+when attempting to fill/interpolate a value at any given location (<em>distance=</em>).
+The distance can be expressed as a number of cells (default) or in the current location's units (if the <em>-m</em> flag is given).
+The latter are typically meters, but can be any other units of a <em>planar</em> coordinate system.
+
+<p>
+
+Note that proper handling of geodetic coordinates (lat/lon) and distances is currently not implemented.
+For lat/lon data, the distance should therefore be specified in cells and usage of <em>r.fill.gaps</em> should
+be restricted to small areas to avoid large inaccuracies that may arise from the different extents of cells
+along the latitudinal and longitudinal axes.
+
+<p>
+
+Distances specfied in map units will be approximated as accurately as the current region's cell resolution
+settings allow. The program will warn if the distance cannot be expressed as whole
+cells at the current region's resolution. In such case, the number of cells in the search window
+will be rounded up. Due to the rounding effect introduced by using cells as spatial units,
+the actual maximum distance considered by the interpolation may be up to half a cell diagonal
+larger than the one specified by the user.
+
+<p>
+
+The interpolator type "wmean" uses a pre-computed matrix of spatial weights To speed up computation. This matrix can
+be examined (printed to the screen) before running the interpolation, by setting the <em>-w</em> flag.
+In mode "wmean", the <em>power=</em> option has the usual meaning: higher values mean
+that cell values in the neighborhood lose their influence on the cell to be interpolated more rapidly with increasing distance from the latter.
+Another way of explaining this effect is to state that larger "power" settings result in more
+localized interpolation, smaller ones in more globalized interpolation. The default setting is <em>power=2.0</em>.
+
+<p>
+
+The interpolators "mean", "median" and "mode" are calculated from all cell values
+within the search radius. No spatial weighting is applied for these methods.
+The "mode" of the input data may not always be unique. In such case, the mode will be the smallest
+value with the highest frequency.
+
+<p>
+
+Often, input data will contain spurious extreme measurements (spikes, outliers, noise) caused by the limits
+of device sensitivity, hardware defects, environmental influences, etc. If the normal, valid
+range of input data is known beforehand, then the <em>minimum=</em> and <em>maximum=</em> options can be used
+to exclude those input cells that have values below or above that range, respectively.
+This will prevent the influence of spikes and outliers from spreading through the interpolation.
+
+<p>
+
+Unless the <em>-p</em> (replace) flag is given, data cells of the input map will be replaced
+with interpolated values instead of preserving them in the output.
+
+Besides the result of the interpolation/gap filling, a second output map can be specified via the
+<em>uncertainty=</em> option. The cell values in this map represent a simple measure of how
+much uncertainty was involved in interpolating each cell value of the primary output map, with
+"0" being the lowest and "1" being the theoretic highest uncertainty. Uncertainty is measured
+by summing up all cells in the neighborhood (defined by the search radius <em>distance=</em>) that contain a value in the <em>input</em> map,
+multiplied by their weights, and dividing the result by the sum of all weights in the neighborhood.
+For <em>mode=wmean</em>, this means that interpolated output cells that
+were computed from many nearby input cells have very low uncertainty and vice versa.
+For all other modes, all weights in the neighborhood are constant "1" and the uncertainty measure
+is a simple measure of how many input data cells were present in the search window.
+
+
+<h2>NOTES</h2>
+
+The key to getting good gap filling results is to understand the spatial weighting scheme
+used in mode "wmean". The weights are precomputed and assigned per cell within the search
+window centered on the location at which to interpolate/gap fill all cells within
+the user-specified distance.
+
+<p>
+
+The illustration below shows the precomputed weights matrix for a search distance of four cells from the center cell:
+
+<pre>
+
+000.00 000.01 000.04 000.07 000.09 000.07 000.04 000.01 000.00
+000.01 000.06 000.13 000.19 000.22 000.19 000.13 000.06 000.01
+000.04 000.13 000.25 000.37 000.42 000.37 000.25 000.13 000.04
+000.07 000.19 000.37 000.56 000.68 000.56 000.37 000.19 000.07
+000.09 000.22 000.42 000.68 001.00 000.68 000.42 000.22 000.09
+000.07 000.19 000.37 000.56 000.68 000.56 000.37 000.19 000.07
+000.04 000.13 000.25 000.37 000.42 000.37 000.25 000.13 000.04
+000.01 000.06 000.13 000.19 000.22 000.19 000.13 000.06 000.01
+000.00 000.01 000.04 000.07 000.09 000.07 000.04 000.01 000.00
+
+</pre>
+
+Note that the weights in such a small window drop rapidly for the default setting of <em>power=2</em>.
+
+<p>
+
+If the distance is given in map units (flag <em>-m</em>), then the search window can be modeled more accurately as a circle.
+The illustration below shows the precomputed weights for a distance in map units that is approximately equivalent to four cells from the center cell:
+
+<pre>
+
+...... ...... ...... 000.00 000.00 000.00 ...... ...... ......
+...... 000.00 000.02 000.06 000.09 000.06 000.02 000.00 ......
+...... 000.02 000.11 000.22 000.28 000.22 000.11 000.02 ......
+000.00 000.07 000.22 000.44 000.58 000.44 000.22 000.07 000.00
+000.00 000.09 000.28 000.58 001.00 000.58 000.28 000.09 000.00
+000.00 000.07 000.22 000.44 000.58 000.44 000.22 000.07 000.00
+...... 000.02 000.11 000.22 000.28 000.22 000.11 000.02 ......
+...... 000.00 000.02 000.06 000.09 000.06 000.02 000.00 ......
+...... ...... ...... 000.00 000.00 000.00 ...... ...... ......
+
+</pre>
+
+<p>
+
+When using a small search radius, <em>cells=</em> must also be set to a small value. Otherwise, there may not be enough cells with data within
+the search radius to support interpolation.
+
+<p>
+
+This module can handle cells with different X and Y resolutions. However, note that the weight matrix
+will be skewed in such cases, with higher weights occurring close to the center and along the axis with the higher resolution.
+This is because weights on the lower resolution axis are less accurately calculated.
+The skewing effect will be stronger if the difference between the X and Y axis resolution is greater and a larger "power" setting is chosen.
+This property of the weights matrix directly reflects the higher information density along the higher resolution axis.
+
+<p>
+
+Note on printing the weights matrix (using the <em>-w</em> flag): the matrix cannot be printed if it is too large.
+
+<p>
+
+The memory estimate provided by the <em>-m</em> flag is a lower limit on the amount of RAM that will be needed.
+
+<p>
+
+If the <em>-s</em> flag is set, floating point type output will be saved as a "single precision" raster map, saving ~50% disk space compared to the default "double precision" output.
+
+
+<h2>EXAMPLES</h2>
+
+Gap-fill a dataset using spatially weighted mean (IDW) and a maximum search radius of 3.0 map units;
+also produce uncertainty estimation map:
+
+<pre>
+r.fill.gaps input=measurements output=result dist=3.0 -m mode=wmean uncertainty=uncert_map
+</pre>
+
+Run a fast low-pass filter (replacement all cells with mean value of neighboring cells) on the input map:
+<pre>
+r.fill.gaps input=measurements output=result dist=10 mode=mean
+</pre>
+
+Fill data gaps in a categorized raster map; preserve existing data:
+<pre>
+r.fill.gaps input=categories output=result dist=100 -m mode=mode -p
+</pre>
+
+
+<h2>SEE ALSO</h2>
+
+<em><a href="http://en.wikipedia.org/wiki/Inverse_distance_weighting">Wikipedia on Inverse Distance Weighting</a></em>,<br>
+<em><a href="r.neighbors.html">r.neighbors</a></em>,<br>
+<em><a href="r.surf.idw2.html">r.surf.idw</a></em>,<br>
+<em><a href="v.surf.idw.html">v.surf.idw</a></em>,<br>
+<em><a href="v.surf.rst.html">v.surf.rst</a></em>
+
+
+<h2>CAVEATS</h2>
+
+The straight-line metric used for converting distances in map units to cell numbers is only adequate for planar coordinate systems.
+Using this module with lat/lon input data will likely give inaccurate results, especially when interpolating over large geographical areas.
+
+<p>
+
+If the distance is set to a relatively large value, processing time
+will quickly approach and eventually exceed that of point-based interpolation modules such
+as <em><a href="v.surf.rst.html">v.surf.rst</a></em>.
+
+
+<h2>AUTHOR</h2>
+
+Benjamin Ducke <br>
+
+<p>
+Latest revision: <i>Mar 20 2015</i>.
Added: grass-addons/grass6/raster/r.fill.gaps/main.c
===================================================================
--- grass-addons/grass6/raster/r.fill.gaps/main.c (rev 0)
+++ grass-addons/grass6/raster/r.fill.gaps/main.c 2017-03-21 11:26:32 UTC (rev 70778)
@@ -0,0 +1,1277 @@
+
+/***************************************************************************
+ *
+ * MODULE: r.fill.gaps (commandline)
+ * AUTHOR(S): Benjamin Ducke
+ * PURPOSE: Rapidly fill "no data" cells in a raster map by simple interpolation.
+ *
+ * COPYRIGHT: (C) 2011-2015 by the GRASS Development Team
+ *
+ * This program is free software under the GPL (>=v2)
+ * Read the file COPYING that comes with GRASS for details.
+ ****************************************************************************
+ */
+
+/*
+
+DOCS:
+- docs need a lot of updates!
+- dropped medoid (since we always have an odd number of cells, median=medoid).
+- flag -p(reserve) reversed (preservation is now default)
+- flag for including center cell has been dropped; center always included
+- "method" renamed to "mode"
+- by default, neighborhood size is now assumed to be in cells,
+ unless -m(ap units) is given
+- only -m(ap units) will result in an exact, circular neighborhood;
+ if (default) cells size specification is used, then the neighborhood
+ shape is a less exact rectangle; in that case the neighborhood dimensions
+ will also be different for x and y if the cell dimensions are different
+- lat/lon data is allowed, but distance measure for -m(ap units) is straight line!
+- cell weights are now normalized to be in range [0;1]
+- center cell weight for "wmean" is now "1.0"
+
+BUGS:
+- mode, median and mode produce large areas of "no data" where there is input data!!!
+
+NEXT VERSION:
+- add lat/lon distance and cost map distance measures (lat/lon currently throws an error).
+- allow user to set which cell value should be filled (i.e. interpreted as "no data")
+
+*/
+
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#include <time.h>
+#include <errno.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+
+#include "cell_funcs.h"
+
+
+/* Global variables :-) */
+double **WEIGHTS; /* neighborhood weights */
+double SUM_WEIGHTS; /* sum of weights in all cells (with or without data) of neighborhood */
+
+unsigned long WINDOW_WIDTH = 0; /* dimensions of neighborhood window */
+unsigned long WINDOW_HEIGHT = 0;
+unsigned long DATA_WIDTH = 0;
+unsigned long DATA_HEIGHT = 0;
+unsigned long PADDING_WIDTH = 0;
+unsigned long PADDING_HEIGHT = 0;
+
+void **CELL_INPUT = NULL;
+void **CELL_INPUT_HANDLES = NULL;
+void *CELL_OUTPUT = NULL;
+FCELL *ERR_OUTPUT = NULL;
+
+/* holds statistics of cells within the neighborhood */
+typedef struct
+{
+ unsigned long num_values; /* number of cells with values in input raster */
+ double *values; /* individual values of all cells */
+ double *weights; /* individual weights of all cells */
+ double result; /* weighted mean of values in entire neighborhood */
+ double certainty; /* certainty measure, always between 0 (lowest) and 1 (highest) */
+ unsigned long *frequencies; /* frequency count for each value */
+ double *overwrite; /* will be set to non-null to overwrite the statistical result with this value */
+} stats_struct;
+
+
+/* function pointers for operation modes */
+void (*GET_STATS) (unsigned long, unsigned long, double, double, int, stats_struct*);
+void (*COLLECT_DATA) (double, double, double, double, stats_struct*);
+
+
+/*
+ * Returns a rough estimate of the amount of RAM
+ * (in bytes) that this program will require.
+ */
+//TODO: this also needs to take into account additional input and output maps (once implemented).
+long int estimate_mem_needed ( long int cols, char *mode )
+{
+ long int mem_count=0;
+ long int in_bytes=0;
+ long int out_bytes=0;
+ long int stat_bytes=0;
+
+
+ /* memory for neighborhood weights and statistics */
+ if ( !strcmp (mode, "wmean" )) {
+ stat_bytes += sizeof (double) * (DATA_WIDTH*DATA_HEIGHT); /* weights matrix */
+ }
+ stat_bytes += sizeof (double) * (DATA_WIDTH*DATA_HEIGHT); /* max. cell values */
+ stat_bytes += sizeof (double) * (DATA_WIDTH*DATA_HEIGHT); /* max. cell weights */
+ stat_bytes += sizeof (unsigned long) * (DATA_WIDTH*DATA_HEIGHT); /* max. cell frequencies */
+
+ /* input data rows with padded buffers */
+ in_bytes = (unsigned long) (WINDOW_HEIGHT * (cols+(PADDING_WIDTH*2)) );
+ if ( IN_TYPE == CELL_TYPE ) {
+ in_bytes += in_bytes * sizeof (CELL);
+ }
+ if ( IN_TYPE == FCELL_TYPE ) {
+ in_bytes += in_bytes * sizeof (FCELL);
+ }
+ if ( IN_TYPE == DCELL_TYPE ) {
+ in_bytes += in_bytes * sizeof (DCELL);
+ }
+
+ /* output data row */
+ out_bytes = (unsigned long) cols;
+ if ( OUT_TYPE == CELL_TYPE ) {
+ out_bytes += out_bytes * sizeof (CELL);
+ }
+ if ( OUT_TYPE == FCELL_TYPE ) {
+ out_bytes += out_bytes * sizeof (FCELL);
+ }
+ if ( OUT_TYPE == DCELL_TYPE ) {
+ out_bytes += out_bytes * sizeof (DCELL);
+ }
+
+ mem_count = stat_bytes + in_bytes + out_bytes;
+
+ return (mem_count);
+}
+
+
+/*
+ * Prints the spatial weights matrix to the console.
+ * This uses a fixed layout which may not be able to print very
+ * large matrices.
+ */
+void print_weights_matrix ( long int rows, long int cols )
+{
+ int i,j;
+ int weight_matrix_line_length = 80;
+ int weight_matrix_weight_length = 7;
+ char weight_matrix_line_buf[weight_matrix_line_length+1];
+ char weight_matrix_weight_buf[weight_matrix_line_length+1];
+
+ G_message (_("Spatial weights neighborhood (cells):"));
+ for ( i = 0; i < rows; i ++ ) {
+ G_strncpy (weight_matrix_line_buf, "", 1);
+ for ( j = 0; j < cols; j ++ ) {
+ if ( WEIGHTS[i][j] != -1.0 ) {
+ snprintf (weight_matrix_weight_buf, weight_matrix_line_length, "%06.2f ", WEIGHTS[i][j]);
+ if ( strlen (weight_matrix_weight_buf) > (weight_matrix_weight_length) ) {
+ snprintf (weight_matrix_weight_buf, weight_matrix_line_length, "[????] ", WEIGHTS[i][j]);
+ }
+ } else {
+ snprintf (weight_matrix_weight_buf, weight_matrix_line_length, "...... ");
+ }
+ if ( strlen (weight_matrix_weight_buf) + strlen (weight_matrix_line_buf) > weight_matrix_line_length ) {
+ G_strncpy ( weight_matrix_line_buf, "[line too long to print]", weight_matrix_line_length);
+ break;
+ } else {
+ G_strcat (weight_matrix_line_buf, weight_matrix_weight_buf);
+ }
+ }
+ fprintf (stdout, "%s\n", weight_matrix_line_buf);
+ }
+}
+
+
+/*
+ * Returns a void* that points to the first valid data cell in the
+ * input data row corresponding to "row_idx".
+ */
+void *get_input_row ( unsigned long row_idx )
+{
+ unsigned long i;
+ void *my_cell=NULL;
+
+
+ my_cell = CELL_INPUT_HANDLES[row_idx];
+
+ for ( i=0; i < PADDING_WIDTH; i ++ )
+ my_cell += CELL_IN_SIZE;
+
+
+ return ( my_cell );
+}
+
+
+/* NEIGHBORHOOD STATISTICS
+ *
+ * The following function provide the different neighborhood statistics (interpolators)
+ * that have been implemented in this software.
+ *
+ * Only those cells go into the statistics that are within the circular neighborhood window.
+ *
+ * Since the input buffers are padded with NULL data for rows lying to the South or North,
+ * West or East of the current region, these functions can just blindly read values above and
+ * below the current position on the W-E axis.
+ *
+ * The parameter "col" must be set to actual GRASS region column number on the W-E axis.
+ *
+ * "min" and "max" define the smallest and largest values to use for interpolation. Set
+ * filter_min and filter_max to !=0 to use these for range filtering the data.
+ *
+ * The results will be stored in the cell_stats object passed to this function. This object
+ * must have been properly initialized before passing it to any of the functions below!
+ */
+
+/*
+ * The different types of neighborhood statistics required different
+ * types of information to be collected.
+ */
+
+void collect_values_unfiltered ( double val1, double val2, double min, double max, stats_struct *stats )
+{
+ stats->values[stats->num_values]=val1;
+ stats->certainty+=val2;
+ stats->num_values++;
+}
+
+void collect_values_filtered ( double val1, double val2, double min, double max, stats_struct *stats)
+{
+ unsigned long i;
+
+ if ( val1 >= min && val1 <= max ) {
+ collect_values_unfiltered (val1,val2,min,max,stats);
+ }
+}
+
+void collect_values_and_weights_unfiltered ( double val1, double val2, double min, double max, stats_struct *stats )
+{
+ stats->values[stats->num_values]=val1;
+ stats->weights[stats->num_values]=val2;
+ stats->certainty+=val2;
+ stats->num_values++;
+}
+
+void collect_values_and_weights_filtered ( double val1, double val2, double min, double max, stats_struct *stats )
+{
+ unsigned long i;
+
+ if ( val1 >= min && val1 <= max ) {
+ collect_values_and_weights_unfiltered (val1,val2,min,max,stats);
+ }
+}
+
+void collect_values_and_frequencies_unfiltered ( double val1, double val2, double min, double max, stats_struct *stats )
+{
+ unsigned long i;
+ int done;
+
+ stats->certainty+=val2;
+
+ /* extreme case: no values collected yet */
+ if ( stats->num_values == 0 ) {
+ stats->values[0]=val1;
+ stats->frequencies[0]=1;
+ stats->num_values++;
+ return;
+ }
+
+ /* at least one value already collected */
+ for ( i=0; i < stats->num_values ; i++ ) {
+ /* does this value already exist in the stats object? */
+ if ( stats->values[i]==val1 ) {
+ /* yes: increase its counter and abort search */
+ stats->frequencies[i]++;
+ stats->values[stats->num_values]=val1;
+ stats->num_values++;
+ return;
+ }
+ }
+ /* no: first occurrence of this value: store as new entry */
+ stats->values[i]=val1;
+ stats->frequencies[i]=1;
+ stats->num_values++;
+}
+
+void collect_values_and_frequencies_filtered ( double val1, double val2, double min, double max, stats_struct *stats )
+{
+ unsigned long i;
+
+ if ( val1 >= min && val1 <= max ) {
+ collect_values_and_frequencies_unfiltered (val1,val2,min,max,stats);
+ }
+}
+
+
+/*
+ * Simple double comparision function for use by qsort().
+ * This is needed for calculating median statistics.
+ */
+int compare_dbl (const void* val1, const void* val2) {
+ if ( *(double*) val1 == *(double*) val2 ) return 0;
+ if ( *(double*) val1 < *(double*) val2 ) return -1;
+ return 1;
+}
+
+/*
+ * Collecting the cell data from the neighborhood is the same
+ * basic loop for every type of statistics: collect only non-null
+ * cells, and collect only within the neighborhood mask.
+ */
+void read_neighborhood (unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct *stats )
+{
+ unsigned long i,j;
+ void *cell;
+ double cell_value;
+
+ /* read data */
+ unsigned long row_position = row_index-PADDING_HEIGHT;
+ stats->num_values=0;
+ stats->certainty=0.0;
+ for ( i = 0; i < DATA_HEIGHT; i ++ ) {
+ cell = CELL_INPUT_HANDLES[i+row_position];
+ cell += CELL_IN_SIZE * col;
+
+ for ( j = 0; j < DATA_WIDTH; j ++ ) {
+ /* read cell from input buffer */
+ if ( !IS_NULL(cell) ) {
+ if ( preserve == TRUE && i == (DATA_HEIGHT-1)/2 && j == (DATA_WIDTH-1)/2 ) {
+ /* center cell: copy through if needed */
+ stats->overwrite = G_malloc (sizeof (double));
+ *stats->overwrite = (double)G_get_raster_value_d(cell,IN_TYPE);
+ } else {
+ /* only add non-null cells to stats */
+ cell_value = (double)G_get_raster_value_d(cell,IN_TYPE);
+ /* only add if within neighborhood */
+ if ( WEIGHTS[i][j] != -1.0 ) {
+ /* get data needed for chosen statistic */
+ COLLECT_DATA (cell_value, WEIGHTS[i][j],min,max,stats);
+ }
+ }
+ }
+ /* go to next cell on current row */
+ cell += CELL_IN_SIZE;
+ }
+ }
+}
+
+
+/*
+ * NEIGHBORHOOD STATISTICS FUNCTION WMEAN
+ * Spatially weighted mean.
+ *
+ * The cell values are multiplied by their spatial weights before they are stored.
+ */
+void get_statistics_wmean (unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct *stats)
+{
+ unsigned long i;
+ double total;
+ double total_weight;
+
+ read_neighborhood (row_index, col, min, max, preserve, stats );
+
+ /* compute weighted average of all valid input cells */
+ total=0;
+ total_weight=0;
+ for (i = 0; i < stats->num_values; i ++) {
+ total += stats->values[i]*stats->weights[i];
+ total_weight += stats->weights[i];
+ }
+ stats->result = total / total_weight;
+}
+
+
+/*
+ * NEIGHBORHOOD STATISTICS FUNCTION MEAN
+ * Simple, unweighted mean.
+ */
+void get_statistics_mean (unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct *stats)
+{
+ unsigned long i;
+ double total;
+
+ read_neighborhood (row_index, col, min, max, preserve, stats );
+
+ /* compute total of all valid input cells */
+ total=0;
+ for (i = 0; i < stats->num_values; i ++) {
+ total += stats->values[i];
+ }
+ stats->result = total / ((double) stats->num_values);
+}
+
+
+/*
+ * NEIGHBORHOOD STATISTICS FUNCTION MEDIAN
+ * Simple, unweighted median. For an even number of data points, the median is the
+ * average of the two central elements in the sorted data list.
+ */
+void get_statistics_median (unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct *stats)
+{
+ read_neighborhood (row_index, col, min, max, preserve, stats );
+
+ /* sort list of values */
+ qsort (&stats->values[0], stats->num_values, sizeof (double), &compare_dbl);
+
+ if ( (double)stats->num_values / 2.0 == 0.0 ) {
+ /* even number of elements: result is average of the two central values */
+ stats->result = (stats->values[stats->num_values/2]+stats->values[(stats->num_values/2)+1])/2.0;
+ } else {
+ /* odd number of elements: result is the central element */
+ stats->result = stats->values[(stats->num_values/2)+1];
+ }
+}
+
+
+/*
+ * NEIGHBORHOOD STATISTICS FUNCTION MODE
+ * Simple, unweighted mode. Mathematically, the mode is not always unique. If there is more than
+ * one value with highest frequency, the smallest one is chosen to represent the mode.
+ */
+void get_statistics_mode (unsigned long row_index, unsigned long col,
+ double min, double max, int preserve,
+ stats_struct *stats)
+{
+ unsigned long i;
+ double mode;
+ unsigned long freq;
+
+ read_neighborhood ( row_index, col, min, max, preserve, stats );
+
+ if (stats->num_values < 1)
+ return;
+
+ mode = stats->values[0];
+ freq = stats->frequencies[0];
+ for (i=1; i < stats->num_values; i++) {
+ if ( stats->frequencies[i] > freq ) {
+ mode = stats->values[i];
+ freq = stats->frequencies[i];
+ }
+ }
+ stats->result=mode;
+}
+
+
+/*
+ * Initializes handlers to point to corresponding data rows.
+ */
+void init_handles ()
+{
+ unsigned long i;
+
+ for ( i=0; i < WINDOW_HEIGHT; i ++ ) {
+ CELL_INPUT_HANDLES[i] = CELL_INPUT[i];
+ }
+}
+
+
+/*
+ * Replaces upper-most data row in input buffer with
+ * a new row from disk.
+ * Re-shuffles the data row handlers, so that the new row
+ * becomes the last row in the index and all other rows move
+ * up one. This procedure is a bit complex, but it is a lot faster
+ * to move around a few pointers in memory, than to re-read all
+ * the data rows for the neighborhood every time we go down one
+ * row in the current region.
+ */
+void advance_one_row (int file_desc, long current_row)
+{
+ unsigned long i,j;
+ void *cell_input;
+ static unsigned long replace_row = 0; /* points to the row which will be replaced next */
+ unsigned long replace_pos = 0;
+
+
+ /* the actual replacement position needs to consider the "no data" padding offset, as well */
+ replace_pos = replace_row + PADDING_HEIGHT;
+
+ /* get address of data row to replace */
+ cell_input = CELL_INPUT[replace_pos];
+ for ( i=0; i < PADDING_WIDTH; i ++ )
+ cell_input += CELL_IN_SIZE;
+
+ /* get next row from disk */
+ G_get_raster_row (file_desc, cell_input, current_row+DATA_HEIGHT, IN_TYPE);
+
+ /* re-assign all row handlers below current replacement position */
+ j = PADDING_HEIGHT;
+ for ( i=0; i < DATA_HEIGHT-(replace_row+1); i ++ ) {
+ CELL_INPUT_HANDLES[j] = CELL_INPUT[replace_pos+1+i];
+ j++;
+ }
+
+ /* re-assign all row handlers up to and including replacement position */
+ for ( i=0; i <= replace_row; i ++ ) {
+ CELL_INPUT_HANDLES[j] = CELL_INPUT[PADDING_HEIGHT+i];
+ j++;
+ }
+
+ replace_row ++;
+ if ( replace_row > (DATA_HEIGHT-1) ) {
+ /* start over once end of data area has been reached */
+ replace_row = 0;
+ }
+}
+
+
+/*
+ * Interpolates one row of input data, stores result in CELL_OUTPUT (global var)
+ */
+void interpolate_row ( unsigned long row_index, unsigned long cols,
+ double min, double max, int preserve,
+ unsigned long min_cells,
+ stats_struct *stats, int write_err )
+{
+ unsigned long j;
+ void *cell_input, *cell_output;
+ FCELL *err_output;
+
+ cell_output = CELL_OUTPUT;
+ err_output = ERR_OUTPUT;
+
+ for ( j = 0; j < cols; j++ ) {
+ /* get neighborhood statistics */
+ GET_STATS (row_index, j, min, max, preserve, stats);
+ /* enough reachable cells in input map? */
+ if ( stats->num_values < min_cells ) {
+ SET_NULL ( cell_output, 1 );
+ if ( write_err )
+ G_set_f_null_value (err_output, 1);
+ } else {
+ if ( stats->overwrite != NULL ) {
+ /* write original value into output map */
+ WRITE_DOUBLE_VAL (cell_output, *stats->overwrite);
+ G_free (stats->overwrite);
+ stats->overwrite = NULL;
+ } else {
+ /* write interpolation result into output map */
+ WRITE_DOUBLE_VAL (cell_output, stats->result);
+ }
+ /* write error/uncertainty output map? */
+ if ( write_err ) {
+ G_set_raster_value_f(err_output,(FCELL)1.0-(stats->certainty/SUM_WEIGHTS), FCELL_TYPE);
+ }
+ }
+ /* advance cell pointers by one cell size */
+ cell_output += CELL_OUT_SIZE;
+ err_output ++;
+ }
+}
+
+
+/*
+ * Pre-computes the matrix of spatial weights.
+ * for operation mode "wmean" (spatially weighted mean), "constant" is passed as "0"
+ * and distance-dependent weigths are calculated. For all other modes, "constant" is
+ * passed as "1" and all cells within the circular neighborhood will be set to "1.0".
+ * In both casses, all cells outside the neighborhood will be set to "-1.0".
+ */
+void build_weights_matrix (double radius, double power, double res_x, double res_y, int constant, int use_map_units )
+{
+
+ unsigned long i, j;
+ double p1_x, p1_y, p2_x, p2_y, A, B, C, W;
+ double tolerance;
+
+
+ /* alloc enough mem for weights matrix */
+ WEIGHTS=G_malloc (sizeof (double*) * DATA_HEIGHT);
+ for ( i = 0; i < DATA_HEIGHT; i ++ ) {
+ WEIGHTS[i] = G_malloc (sizeof(double)*DATA_WIDTH);
+ }
+
+ /* center of the neighborhood window in real map units */
+ p1_x = (DATA_WIDTH/2 * res_x) + (res_x / 2.0);
+ p1_y = (DATA_HEIGHT/2 * res_y) + (res_y / 2.0);
+
+ /* tolerance for including half cells in the neighborhood */
+ tolerance = (sqrt (pow(res_x,2)+pow(res_y,2))) / 2;
+
+ /* 1st pass: get largest possible weight for normalization */
+ double max = -1.0;
+ for ( i = 0; i < DATA_HEIGHT; i ++ ) {
+ for ( j = 0; j < DATA_WIDTH; j ++ ) {
+ p2_x = (j * res_x) + (res_x / 2.0);
+ p2_y = (i * res_y) + (res_y / 2.0);
+ A = fabs (p2_x - p1_x);
+ B = fabs (p2_y - p1_y);
+ C = sqrt (pow(A,2)+pow(B,2));
+ if ( use_map_units ) {
+ if (C > radius+tolerance) {
+ WEIGHTS[i][j] = -1.0;
+ } else {
+ WEIGHTS[i][j] = C;
+ }
+ } else {
+ WEIGHTS[i][j] = C;
+ }
+ if ( WEIGHTS[i][j] > max ) {
+ /* update max value */
+ max = WEIGHTS[i][j];
+ }
+ }
+ }
+
+ /* Build the weights matrix */
+ SUM_WEIGHTS = 0.0;
+ for ( i = 0; i < DATA_HEIGHT; i ++ ) {
+ for ( j = 0; j < DATA_WIDTH; j ++ ) {
+ /* Assign neighborhood coordinates with 0/0
+ at the top left of the cell neighborhood matrix. */
+ p2_x = (j * res_x) + (res_x / 2.0);
+ p2_y = (i * res_y) + (res_y / 2.0);
+ /* get distance from window center */
+ A = fabs (p2_x - p1_x);
+ B = fabs (p2_y - p1_y);
+ C = sqrt (pow(A,2)+pow(B,2));
+ if ( constant ) {
+ W = 1.0;
+ } else {
+ W = ((pow (1-(C/max),power)));
+ }
+ /* exclude neighborhood locations that are farther
+ from the center than the user-defined distance
+ plus a tolerance of half the current region's
+ cell diagonal */
+ if ( use_map_units ) {
+ if (C > radius+tolerance) {
+ WEIGHTS[i][j] = -1.0;
+ } else {
+ WEIGHTS[i][j] = W;
+ WEIGHTS[DATA_HEIGHT/2][DATA_WIDTH/2] = 0.0; /* safeguard against total weight growing by center cell weight (InF) */
+ SUM_WEIGHTS += WEIGHTS[i][j];
+ }
+ } else {
+ WEIGHTS[i][j] = W;
+ WEIGHTS[DATA_HEIGHT/2][DATA_WIDTH/2] = 0.0; /* safeguard against total weight growing by center cell weight (InF) */
+ SUM_WEIGHTS += WEIGHTS[i][j];
+ }
+ }
+ }
+
+ /* weight of center cell is always = 1 */
+ WEIGHTS[DATA_HEIGHT/2][DATA_WIDTH/2] = 1.0;
+}
+
+
+/*
+ *
+ * MAIN FUNCTION
+ *
+ */
+int main(int argc, char *argv[])
+{
+ /* processing time keepers */
+ time_t start, finish;
+
+ /* GRASS region properties */
+ struct Cell_head cellhd, region;
+ struct Range int_range;
+ struct FPRange fp_range;
+ struct History hist;
+ unsigned long rows = 0;
+ unsigned long cols = 0;
+ double res_x, res_y = 0.0;
+ int projection;
+
+ /* GRASS module options */
+ struct GModule *module;
+ struct
+ {
+ struct Option
+ *input, *output, *error,
+ *radius, *mode, *power, *min, *max, *minpts;
+ struct Flag
+ *dist_m, *preserve, *print_w, *print_u, *center, *single_precision;
+ } parm;
+
+ /* program settings */
+ char *input;
+ char *output;
+ char *mapset;
+ char *error;
+ double radius=1.0;
+ unsigned long min_cells=12;
+ double power=2.0;
+ double min=0.0;
+ double max=0.0;
+ int filter_min=0;
+ int filter_max=0;
+ int write_error;
+
+ /* file handlers */
+ void *cell_input;
+ void *cell_output;
+ int in_fd;
+ int out_fd;
+ int err_fd;
+
+ /* cell statistics object */
+ stats_struct cell_stats;
+
+ /* generic indices, loop counters, etc. */
+ unsigned long i,j,k;
+ long l;
+ unsigned long count;
+
+
+ start = time(NULL);
+
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->keywords = _("raster,interpolation,IDW");
+ module->description =
+ _("Rapidly fills 'no data' cells of a raster map with interpolated values (IDW).");
+
+ /* parameters */
+
+ parm.input = G_define_standard_option(G_OPT_R_INPUT);
+ parm.input->key = "input";
+ parm.input->required = YES;
+ parm.input->multiple = NO;
+ parm.input->description = _("Raster map with data gaps to fill");
+
+ parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.output->required = YES;
+ parm.output->key = "output";
+ parm.output->description = _("Name of result output map");
+
+ parm.error = G_define_standard_option(G_OPT_R_OUTPUT);
+ parm.error->required = NO;
+ parm.error->key = "uncertainty";
+ parm.error->description = _("Name of uncertainty output map");
+
+ parm.radius = G_define_option();
+ parm.radius->key = "distance";
+ parm.radius->key_desc = "value";
+ parm.radius->required = YES;
+ parm.radius->multiple = NO;
+ parm.radius->type = TYPE_DOUBLE;
+ parm.radius->description = _("Distance threshold (default: in cells) for interpolation");
+ parm.radius->answer = "3";
+
+ parm.mode = G_define_option();
+ parm.mode->key = "mode";
+ parm.mode->key_desc = "name";
+ parm.mode->required = YES;
+ parm.mode->multiple = NO;
+ parm.mode->type = TYPE_STRING;
+ parm.mode->description = _("Statistic for interpolated cell values");
+ parm.mode->options = "wmean,mean,median,mode";
+ parm.mode->answer = "wmean";
+
+ parm.min = G_define_option();
+ parm.min->key = "minimum";
+ parm.min->key_desc = "value";
+ parm.min->required = NO;
+ parm.min->multiple = NO;
+ parm.min->type = TYPE_DOUBLE;
+ parm.min->description = _("Minimum input data value to include in interpolation");
+
+ parm.max = G_define_option();
+ parm.max->key = "maximum";
+ parm.max->key_desc = "value";
+ parm.max->required = NO;
+ parm.max->multiple = NO;
+ parm.max->type = TYPE_DOUBLE;
+ parm.max->description = _("Maximum input data value to include in interpolation");
+
+ parm.power = G_define_option();
+ parm.power->key = "power";
+ parm.power->key_desc = "value";
+ parm.power->required = YES;
+ parm.power->multiple = NO;
+ parm.power->type = TYPE_DOUBLE;
+ parm.power->answer = "2.0";
+ parm.power->description = _("Power coefficient for IDW interpolation");
+
+ parm.minpts = G_define_option();
+ parm.minpts->key = "cells";
+ parm.minpts->key_desc = "value";
+ parm.minpts->required = YES;
+ parm.minpts->multiple = NO;
+ parm.minpts->type = TYPE_INTEGER;
+ parm.minpts->answer = "8";
+ parm.minpts->description = _("Minimum number of data cells within search radius");
+
+ parm.dist_m = G_define_flag();
+ parm.dist_m->key = 'm';
+ parm.dist_m->description = _("Interpret distance as map units, not cell number");
+
+ parm.preserve = G_define_flag();
+ parm.preserve->key = 'p';
+ parm.preserve->description = _("Preserve original cell values");
+
+ parm.print_w = G_define_flag();
+ parm.print_w->key = 'w';
+ parm.print_w->description = _("Just print the spatial weights matrix");
+
+ parm.print_u = G_define_flag();
+ parm.print_u->key = 'u';
+ parm.print_u->description = _("Just print estimated memory usage");
+
+ parm.single_precision = G_define_flag();
+ parm.single_precision->key = 's';
+ parm.single_precision->description = _("Single precision floating point output");
+
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ input = parm.input->answer;
+ output = parm.output->answer;
+
+ /* get setting of current GRASS region */
+ G_get_window(®ion);
+ projection = region.proj;
+ if ( projection == PROJECTION_LL && parm.dist_m->answer ) {
+ G_warning (_("You are working with lat/lon data."));
+ G_warning (_("This module uses a straight-line distance metric."));
+ G_warning (_("Expect inaccuracies."));
+ }
+ rows = (unsigned long) region.rows;
+ cols = (unsigned long) region.cols;
+ res_x = (double) region.ew_res;
+ res_y = (double) region.ns_res;
+
+ /* get user parameters */
+ radius = strtod (parm.radius->answer,0);
+ power = strtod (parm.power->answer,0);
+ min_cells = atol (parm.minpts->answer);
+ write_error = 0;
+ if ( parm.error->answer ) {
+ write_error = 1;
+ if (!strcmp(parm.error->answer,parm.output->answer)) {
+ G_fatal_error (_("Result map name cannot be identical with uncertainty map name."));
+ }
+ }
+
+ /* validate user parameters */
+ double max_dist = 0.0;
+ if ( parm.dist_m->answer ) {
+ if ( radius < 0 ) {
+ G_fatal_error (_("Maximum distance must be larger than zero."));
+ }
+ if ( res_x < res_y ) {
+ if ( radius < res_x )
+ G_fatal_error (_("Maximum distance must be at least '%.6f' (W-E resolution)."), res_x);
+ }
+ if ( res_y < res_x ) {
+ if ( radius < res_y )
+ G_fatal_error (_("Maximum distance must be at least '%.6f' (S-N resolution)."), res_y);
+ }
+ if ( res_y == res_x ) {
+ if ( radius < res_y )
+ G_fatal_error (_("Maximum distance must be at least '%.6f' (W-E and S-N resolution)."), res_y);
+ }
+ max_dist = sqrt ( pow((cols*res_x),2) + pow((rows*res_y),2) );
+ if ( radius > max_dist ) {
+ G_warning (_("Maximum distance too large. Adjusted to '%.6f' (diagonal of current region)."), max_dist);
+ radius = max_dist;
+ }
+ } else {
+ unsigned long radius_i = (unsigned long)radius;
+ radius=(double)radius_i; //truncate to whole cell number
+ if ( radius < 1 ) {
+ G_fatal_error (_("Maximum distance must be at least one cell."));
+ }
+ unsigned long max_dist_i = (unsigned long)(sqrt ( pow((cols),2) + pow((rows),2) ));
+ max_dist = (double) max_dist_i;
+ if ( radius > max_dist ) {
+ G_warning (_("Maximum distance too large. Adjusted to '%lu' cells (diagonal of current region)."), max_dist_i);
+ radius = (double)max_dist_i;
+ }
+ }
+
+ if ( min_cells < 1 ) {
+ G_fatal_error (_("Minimum number of cells must be at least '1'."));
+ }
+ if ( min_cells > ((DATA_WIDTH*DATA_HEIGHT)-1)) {
+ G_fatal_error (_("Specified minimum number of cells unreachable with current settings."));
+ }
+
+ if ( filter_min != 0 && filter_max != 0 ) {
+ if ( min >= max ) {
+ G_fatal_error (_("Value for 'minimum' must be smaller than value for 'maximum'."));
+ }
+ }
+ if ( parm.power->answer && strcmp (parm.mode->answer,"wmean") ) {
+ G_warning (_("The 'power' option has no effect in any mode other than 'wmean'."));
+ parm.power->answer=0;
+ }
+
+
+ /* rounded dimensions of weight matrix (in map cells) */
+ if ( parm.dist_m->answer ) {
+ DATA_WIDTH = ((unsigned long) ceil (radius/res_x)) * 2 + 1;
+ DATA_HEIGHT = ((unsigned long) ceil (radius/res_y)) * 2 + 1;
+ if ((!parm.print_w->answer) && (fmod (radius,res_x) != 0 || fmod (radius,res_y) != 0)) {
+ G_warning (_("The specified maximum distance cannot be resolved to whole cells\n at the current resolution settings."));
+ }
+ } else {
+ DATA_WIDTH = (unsigned long) radius*2 + 1;
+ DATA_HEIGHT = (unsigned long) radius*2 + 1;
+ }
+ PADDING_WIDTH = (DATA_WIDTH-1)/2;
+ PADDING_HEIGHT = (DATA_HEIGHT-1)/2;
+ WINDOW_WIDTH = (PADDING_WIDTH*2) + DATA_WIDTH;
+ WINDOW_HEIGHT = (PADDING_HEIGHT*2) + DATA_HEIGHT;
+ G_message (_("W-E size of neighborhood is %lu cells."), DATA_WIDTH);
+ G_message (_("S-N size of neighborhood is %lu cells."), DATA_HEIGHT);
+
+ if ( DATA_WIDTH < 3 || DATA_HEIGHT < 3 ) {
+ G_fatal_error (_("Neighborhood cannot be smaller than 3 cells in X or Y direction."));
+ }
+
+ if (parm.print_w->answer) {
+ if ( !strcmp (parm.mode->answer,"wmean") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 0, parm.dist_m->answer);
+ }
+ if ( !strcmp (parm.mode->answer,"mean") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ }
+ if ( !strcmp (parm.mode->answer,"median") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ }
+ if ( !strcmp (parm.mode->answer,"mode") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ }
+ print_weights_matrix (DATA_HEIGHT, DATA_WIDTH);
+ /* continue only if "-u" flag has also been given */
+ if ( !parm.print_u->answer )
+ exit (0);
+ }
+
+ /* open raster input map and get its storage type */
+ mapset = G_find_cell2(input, "");
+ if (!mapset)
+ G_fatal_error(_("Raster map <%s> not found"), input);
+ if (G_get_cellhd(input, mapset, &cellhd) < 0)
+ G_fatal_error(_("Unable to read header of raster map <%s@%s>"), input, mapset);
+ if ((in_fd = G_open_cell_old(input, mapset)) < 0)
+ G_fatal_error(_("Unable to open raster map <%s> in mapset <%s>"), input, mapset);
+ IN_TYPE = G_get_raster_map_type(in_fd);
+
+ /* minimum and maximum values for interpolating range */
+ if ( IN_TYPE == CELL_TYPE ) {
+ G_read_range ( input, mapset, &int_range);
+ min = (double) int_range.min;
+ max = (double) int_range.max;
+ } else {
+ G_read_fp_range ( input, mapset, &fp_range);
+ min = (double) fp_range.min;
+ max = (double) fp_range.max;
+ }
+ if ( parm.min->answer ) {
+ min = strtod (parm.min->answer,0);
+ filter_min = 1;
+ }
+ if ( parm.max->answer ) {
+ max = strtod (parm.max->answer,0);
+ filter_max = 1;
+ }
+ G_message (_("Input data range is %f to %f.\n"), min, max);
+
+ /* determine input and output data types, advise user */
+ OUT_TYPE = IN_TYPE;
+ if ( IN_TYPE == DCELL_TYPE ) {
+ if ( parm.single_precision->answer ) {
+ OUT_TYPE = FCELL_TYPE;
+ } else {
+ OUT_TYPE = DCELL_TYPE;
+ }
+ }
+ if ( IN_TYPE == CELL_TYPE ) {
+ if ( (!strcmp (parm.mode->answer,"wmean")) || (!strcmp (parm.mode->answer,"mean"))
+ || (!strcmp (parm.mode->answer,"median"))) {
+ G_warning (_("Input data type is integer but interpolation mode is '%s'."), parm.mode->answer);
+ if ( parm.single_precision->answer ) {
+ OUT_TYPE = FCELL_TYPE;
+ G_warning (_("Output type changed to floating point (single)."));
+ } else {
+ OUT_TYPE = DCELL_TYPE;
+ G_warning (_("Output type changed to floating point (double)."));
+ }
+ } else {
+ if ( parm.single_precision->answer ) {
+ G_warning (_("Ignoring '%c' flag. Output data type will be integer."), parm.single_precision->key);
+ }
+ }
+ }
+ char data_type_string_in[21];
+ char data_type_string_out[21];
+ if ( IN_TYPE == CELL_TYPE ) {
+ G_strncpy (data_type_string_in, "integer", 20);
+ }
+ else if ( IN_TYPE == FCELL_TYPE ) {
+ G_strncpy (data_type_string_in, "single", 20);
+ }
+ else if ( IN_TYPE == DCELL_TYPE ) {
+ G_strncpy (data_type_string_in, "double", 20);
+ }
+ if ( OUT_TYPE == CELL_TYPE ) {
+ G_strncpy (data_type_string_out, "integer", 20);
+ }
+ else if ( OUT_TYPE == FCELL_TYPE ) {
+ G_strncpy (data_type_string_out, "single", 20);
+ }
+ else if ( OUT_TYPE == DCELL_TYPE ) {
+ G_strncpy (data_type_string_out, "double", 20);
+ }
+
+ /* initialize data type dependent cell handling functions */
+ init_cell_funcs();
+
+ G_message (_("Input data type is '%s' (%i bytes) and output data type is '%s' (%i bytes)."),
+ data_type_string_in, CELL_IN_SIZE, data_type_string_out, CELL_OUT_SIZE );
+
+ /* just print projected mem usage if user wants it so */
+ G_message ("Minimal estimated memory usage is %.3f MB.", ((double) estimate_mem_needed(cols,parm.mode->answer))/1024/1024);
+ if ( parm.print_u->answer ) {
+ exit (0);
+ }
+
+ /* Allocate enough memory to read n="distance"x2+1 rows of input map data,
+ * plus a buffer of "distance" size above and below the actual data
+ * rows, and to the right and left.
+ * The buffer will be filled with "no data" cells and has the
+ * effect that we can later read data anywhere within the search
+ * neighborhood, without having to worry about alignment problems.
+ * */
+ CELL_INPUT = G_malloc ( CELL_IN_PTR_SIZE * WINDOW_HEIGHT );
+ for ( i = 0; i < WINDOW_HEIGHT; i ++ ) {
+ CELL_INPUT[i]= G_malloc ( CELL_IN_SIZE * (cols + (PADDING_WIDTH*2)));
+ }
+ for ( i = 0; i < WINDOW_HEIGHT; i ++ ) {
+ G_set_null_value (CELL_INPUT[i], cols + (PADDING_WIDTH*2), IN_TYPE );
+ }
+
+ /*
+ * Allocate array of raster row data handlers.
+ * When reading rows from the input data buffer, we can use these
+ * handlers instead of the original data pointers. That way, we
+ * only need to read one new row of data from the disk each time
+ * the neighborhood advances down one row in the region. Then, we
+ * re-shuffle the row handlers, so that the first handler always
+ * points to the first row of the neighborhood data, and the following
+ * ones to the subsequent rows.
+ * This should save a lot of disk access time.
+ */
+ CELL_INPUT_HANDLES = G_malloc (CELL_IN_PTR_SIZE*WINDOW_HEIGHT);
+
+ /* create statistics object */
+ cell_stats.values=G_malloc (sizeof(double)*WINDOW_WIDTH*WINDOW_HEIGHT);
+ cell_stats.weights=G_malloc (sizeof(double)*WINDOW_WIDTH*WINDOW_HEIGHT);
+ cell_stats.frequencies=G_malloc (sizeof(unsigned long)*WINDOW_WIDTH*WINDOW_HEIGHT);
+ cell_stats.overwrite = NULL;
+
+ /* set statistics functions according to user option setting */
+ if ( !strcmp (parm.mode->answer,"wmean") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 0, parm.dist_m->answer);
+ GET_STATS = &get_statistics_wmean;
+ if ( filter_min == 1 || filter_max == 1 ) {
+ COLLECT_DATA=&collect_values_and_weights_filtered;
+ } else {
+ COLLECT_DATA=&collect_values_and_weights_unfiltered;
+ }
+ }
+ if ( !strcmp (parm.mode->answer,"mean") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ GET_STATS = &get_statistics_mean;
+ if ( filter_min == 1 || filter_max == 1 ) {
+ COLLECT_DATA=&collect_values_filtered;
+ } else {
+ COLLECT_DATA=&collect_values_unfiltered;
+ }
+ }
+ if ( !strcmp (parm.mode->answer,"median") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ GET_STATS = &get_statistics_median;
+ if ( filter_min == 1 || filter_max == 1 ) {
+ COLLECT_DATA=&collect_values_and_frequencies_filtered;
+ } else {
+ COLLECT_DATA=&collect_values_and_frequencies_unfiltered;
+ }
+ }
+ if ( !strcmp (parm.mode->answer,"mode") ) {
+ build_weights_matrix (radius, power, res_x, res_y, 1, parm.dist_m->answer);
+ GET_STATS = &get_statistics_mode;
+ if ( filter_min == 1 || filter_max == 1 ) {
+ COLLECT_DATA=&collect_values_and_frequencies_filtered;
+ } else {
+ COLLECT_DATA=&collect_values_and_frequencies_unfiltered;
+ }
+ }
+
+ /*
+ *
+ * MAIN LOOP
+ *
+ */
+
+ /* Open output map with right data type */
+ out_fd = G_open_raster_new(output, OUT_TYPE);
+ if (out_fd < 0) {
+ G_fatal_error ("Cannot open output map.");
+ exit(EXIT_FAILURE);
+ }
+
+ /* Reserve memory for one output row buffer */
+ CELL_OUTPUT=G_allocate_raster_buf (OUT_TYPE);
+
+ /* initialize output row */
+ SET_NULL (CELL_OUTPUT, cols);
+
+ /* produce uncertainty output map? */
+ if ( parm.error->answer ) {
+ /* Open output map with right data type */
+ err_fd = G_open_raster_new(parm.error->answer, FCELL_TYPE);
+ if (err_fd < 0) {
+ G_fatal_error ("Cannot open uncertainty output map.");
+ exit(EXIT_FAILURE);
+ }
+ ERR_OUTPUT=G_allocate_raster_buf (FCELL_TYPE);
+ /* initialize output row */
+ G_set_f_null_value (ERR_OUTPUT, cols);
+ }
+
+ /* row indices to handle input data buffer */
+ unsigned long center_row = (PADDING_HEIGHT*2);
+ unsigned long row_idx=0;
+
+ /* Visit every row in the input dataset.
+ * To avoid making this code complex and having a lot of
+ * if/then boundary checks while looping, processing has
+ * been split into upper edge, main part and lower edge.
+ * The upper and lower edge comprise the diameter of the
+ * neighborhood window.
+ * */
+
+ G_message (_("Interpolating:"));
+ unsigned long current_row=0;
+
+ /* first part: upper edge of region */
+ init_handles();
+ for ( i=0; i < DATA_HEIGHT; i++ ) {
+ cell_input = get_input_row(PADDING_HEIGHT+i);
+
+ cell_input = CELL_INPUT[PADDING_HEIGHT+i];
+ for (j=0; j < PADDING_WIDTH; j ++) {
+ cell_input += CELL_IN_SIZE;
+ }
+ G_get_raster_row (in_fd, cell_input, i, IN_TYPE);
+ }
+ for ( i = 0; i <= PADDING_HEIGHT; i++ ) {
+ row_idx = PADDING_HEIGHT + i;
+ interpolate_row ( row_idx, cols, min, max, parm.preserve->answer, min_cells,
+ &cell_stats, write_error );
+ /* write output row buffer to disk */
+ G_put_raster_row (out_fd, CELL_OUTPUT, OUT_TYPE);
+ if ( parm.error->answer )
+ G_put_raster_row (err_fd, ERR_OUTPUT, FCELL_TYPE);
+ G_percent(current_row+1, rows, 2);
+ current_row++;
+ }
+
+ /* second part: region between upper and lower edge */
+ for ( i = 0; i < rows - (DATA_HEIGHT+1); i++ ) {
+ l = i;
+ advance_one_row (in_fd,l);
+ l++;
+ row_idx = center_row;
+ interpolate_row ( row_idx, cols, min, max, parm.preserve->answer, min_cells,
+ &cell_stats, write_error );
+ /* write output row buffer to disk */
+ G_put_raster_row (out_fd, CELL_OUTPUT, OUT_TYPE);
+ if ( parm.error->answer )
+ G_put_raster_row (err_fd, ERR_OUTPUT, FCELL_TYPE);
+ G_percent(current_row+1, rows, 2);
+ current_row++;
+ }
+
+ /* third part: lower edge */
+ init_handles();
+ for ( i=rows - DATA_HEIGHT ; i < rows; i++ ) {
+ row_idx = (DATA_HEIGHT+PADDING_HEIGHT) - (rows-i);
+ cell_input = get_input_row(row_idx);
+ G_get_raster_row (in_fd, cell_input, i, IN_TYPE);
+ }
+ for ( i = rows - PADDING_HEIGHT-1; i < rows; i++ ) {
+ row_idx = PADDING_HEIGHT + (DATA_HEIGHT) - (rows-i);
+ interpolate_row (row_idx, cols, min, max, parm.preserve->answer, min_cells,
+ &cell_stats, write_error );
+ /* write output row buffer to disk */
+ G_put_raster_row (out_fd, CELL_OUTPUT, OUT_TYPE);
+ if ( parm.error->answer )
+ G_put_raster_row (err_fd, ERR_OUTPUT, FCELL_TYPE);
+ G_percent(current_row+1, rows, 2);
+ current_row++;
+ }
+
+ /* close all maps */
+ G_close_cell(out_fd);
+ G_close_cell(in_fd);
+ if ( parm.error->answer ) {
+ G_close_cell(err_fd);
+ }
+
+ /* Free memory */
+ for ( i = 0; i < DATA_HEIGHT; i++ ) {
+ G_free (WEIGHTS[i]);
+ }
+ G_free (WEIGHTS);
+
+ if ( CELL_INPUT != NULL ) {
+ for ( i = 0; i < DATA_HEIGHT; i ++ ) {
+ G_free (CELL_INPUT[i]);
+ }
+ G_free (CELL_INPUT);
+ }
+
+ if ( CELL_OUTPUT != NULL )
+ G_free (CELL_OUTPUT);
+
+ if ( parm.error->answer )
+ G_free (ERR_OUTPUT);
+
+ G_free (cell_stats.values);
+ G_free (cell_stats.weights);
+ G_free (cell_stats.frequencies);
+
+ /* write metadata into result and error maps */
+ G_short_history (parm.output->answer,"raster",&hist);
+ sprintf (hist.title,"Result of interpolation/gap filling");
+ hist.edlinecnt=2;
+ if ( parm.dist_m->answer ) {
+ sprintf (hist.edhist[0],"Settings: mode=%s, distance (map units)=%.6f, power=%.3f",
+ parm.mode->answer, radius, power );
+ } else {
+ sprintf (hist.edhist[0],"Settings: mode=%s, distance (cells)=%lu, power=%.3f",
+ parm.mode->answer, (unsigned long) radius, power );
+ }
+ sprintf (hist.edhist[1]," min=%.3f, max=%.3f, min. points=%lu",
+ min, max, min_cells );
+ G_write_history (parm.output->answer, &hist);
+
+ if ( parm.error->answer ) {
+ G_short_history (parm.error->answer,"raster",&hist);
+ sprintf (hist.title,"Uncertainty of interpolation/gap filling");
+ hist.edlinecnt=2;
+ sprintf (hist.edhist[0],"Result map: %s", parm.output->answer );
+ sprintf (hist.edhist[1],"Theoretic range is '0' (lowest) to '1' (highest).");
+ G_write_history (parm.error->answer, &hist);
+ }
+
+ finish = time(NULL);
+ double ticks = difftime(finish,start);
+ int hours = trunc (ticks/3600);
+ ticks = ticks - (hours*3600);
+ int minutes = trunc (ticks/60);
+ ticks = ticks - (minutes*60);
+ int seconds = trunc (ticks);
+
+ /* DONE */
+ G_done_msg(_("Processing time was %ih%im%is."),hours,minutes,seconds);
+
+ return (EXIT_SUCCESS);
+}
Added: grass-addons/grass6/raster/r.fill.gaps/r.fill.gaps.01.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass6/raster/r.fill.gaps/r.fill.gaps.01.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass6/raster/r.fill.gaps/r.fill.gaps.02.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass6/raster/r.fill.gaps/r.fill.gaps.02.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
More information about the grass-commit
mailing list