[GRASS-SVN] r70880 - grass-addons/grass7/raster/r.fill.gaps
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Apr 14 15:37:30 PDT 2017
Author: wenzeslaus
Date: 2017-04-14 15:37:29 -0700 (Fri, 14 Apr 2017)
New Revision: 70880
Added:
grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.html
grass-addons/grass7/raster/r.fill.gaps/r_fill_gaps_01.png
grass-addons/grass7/raster/r.fill.gaps/r_fill_gaps_02.png
Removed:
grass-addons/grass7/raster/r.fill.gaps/description.html
grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.01.png
grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.02.png
Modified:
grass-addons/grass7/raster/r.fill.gaps/Makefile
grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c
grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h
grass-addons/grass7/raster/r.fill.gaps/main.c
Log:
r.fill.gaps: port code from G6 to G7
Modified: grass-addons/grass7/raster/r.fill.gaps/Makefile
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/Makefile 2017-04-14 22:31:39 UTC (rev 70879)
+++ grass-addons/grass7/raster/r.fill.gaps/Makefile 2017-04-14 22:37:29 UTC (rev 70880)
@@ -2,10 +2,9 @@
PGM = r.fill.gaps
-LIBES = $(GISLIB)
-DEPENDENCIES = $(GISDEP)
+LIBES = $(GISLIB) $(RASTERLIB)
+DEPENDENCIES = $(GISDEP) $(RASTERDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
default: cmd
-
Modified: grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c 2017-04-14 22:31:39 UTC (rev 70879)
+++ grass-addons/grass7/raster/r.fill.gaps/cell_funcs.c 2017-04-14 22:37:29 UTC (rev 70880)
@@ -21,31 +21,46 @@
#include <math.h>
#include <string.h>
#include <stdlib.h>
+
#include <grass/gis.h>
#include <grass/glocale.h>
#include "cell_funcs.h"
+RASTER_MAP_TYPE IN_TYPE;
+RASTER_MAP_TYPE OUT_TYPE;
+
+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 *);
+void (*WRITE_DOUBLE_VAL) (void *, double);
+int (*IS_NULL) (void *);
+void (*SET_NULL) (void *, unsigned long);
+
/*
* 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);
+ Rast_set_c_value(cell_output, Rast_get_c_value(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);
+ Rast_set_f_value(cell_output, Rast_get_f_value(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);
+ Rast_set_d_value(cell_output, Rast_get_d_value(cell_input, IN_TYPE),
+ OUT_TYPE);
}
@@ -55,17 +70,17 @@
void write_double_value_c(void *cell, double val)
{
- G_set_raster_value_c(cell, (CELL) val, OUT_TYPE);
+ Rast_set_c_value(cell, (CELL) val, OUT_TYPE);
}
void write_double_value_f(void *cell, double val)
{
- G_set_raster_value_f(cell, (FCELL) val, OUT_TYPE);
+ Rast_set_f_value(cell, (FCELL) val, OUT_TYPE);
}
void write_double_value_d(void *cell, double val)
{
- G_set_raster_value_d(cell, (DCELL) val, OUT_TYPE);
+ Rast_set_d_value(cell, (DCELL) val, OUT_TYPE);
}
@@ -75,17 +90,17 @@
*/
int is_null_value_c(void *cell)
{
- return (G_is_c_null_value((CELL *) cell));
+ return (Rast_is_c_null_value((CELL *) cell));
}
int is_null_value_f(void *cell)
{
- return (G_is_f_null_value((FCELL *) cell));
+ return (Rast_is_f_null_value((FCELL *) cell));
}
int is_null_value_d(void *cell)
{
- return (G_is_d_null_value((DCELL *) cell));
+ return (Rast_is_d_null_value((DCELL *) cell));
}
@@ -94,17 +109,17 @@
*/
void set_null_c(void *cells, unsigned long count)
{
- G_set_c_null_value((CELL *) cells, count);
+ Rast_set_c_null_value((CELL *) cells, count);
}
void set_null_f(void *cells, unsigned long count)
{
- G_set_f_null_value((FCELL *) cells, count);
+ Rast_set_f_null_value((FCELL *) cells, count);
}
void set_null_d(void *cells, unsigned long count)
{
- G_set_d_null_value((DCELL *) cells, count);
+ Rast_set_d_null_value((DCELL *) cells, count);
}
Modified: grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h 2017-04-14 22:31:39 UTC (rev 70879)
+++ grass-addons/grass7/raster/r.fill.gaps/cell_funcs.h 2017-04-14 22:37:29 UTC (rev 70880)
@@ -17,19 +17,23 @@
#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;
+#include <grass/raster.h>
+void init_cell_funcs();
+
+extern RASTER_MAP_TYPE IN_TYPE; /* stuff for cell input and output data */
+extern 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;
+extern unsigned char CELL_IN_SIZE;
+extern unsigned char CELL_IN_PTR_SIZE;
+extern unsigned char CELL_OUT_SIZE;
+extern unsigned char CELL_OUT_PTR_SIZE;
+extern 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) */
+extern void (*WRITE_CELL_VAL) (void *, void *); /* write a cell value of any type into an output cell */
+extern void (*WRITE_DOUBLE_VAL) (void *, double); /* write a double value into an output cell */
+extern int (*IS_NULL) (void *); /* check if a cell is "no data" */
+extern void (*SET_NULL) (void *, unsigned long); /* write null value(s) */
#endif /* CELL_FUNCS_H */
Deleted: grass-addons/grass7/raster/r.fill.gaps/description.html
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/description.html 2017-04-14 22:31:39 UTC (rev 70879)
+++ grass-addons/grass7/raster/r.fill.gaps/description.html 2017-04-14 22:37:29 UTC (rev 70880)
@@ -1,255 +0,0 @@
-<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>.
Modified: grass-addons/grass7/raster/r.fill.gaps/main.c
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/main.c 2017-04-14 22:31:39 UTC (rev 70879)
+++ grass-addons/grass7/raster/r.fill.gaps/main.c 2017-04-14 22:37:29 UTC (rev 70880)
@@ -152,7 +152,7 @@
G_message(_("Spatial weights neighborhood (cells):"));
for (i = 0; i < rows; i++) {
- G_strncpy(weight_matrix_line_buf, "", 1);
+ weight_matrix_line_buf[0] = '\0';
for (j = 0; j < cols; j++) {
if (WEIGHTS[i][j] != -1.0) {
snprintf(weight_matrix_weight_buf, weight_matrix_line_length,
@@ -170,12 +170,12 @@
}
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);
+ 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);
+ strcat(weight_matrix_line_buf, weight_matrix_weight_buf);
}
}
fprintf(stdout, "%s\n", weight_matrix_line_buf);
@@ -356,11 +356,11 @@
/* center cell: copy through if needed */
stats->overwrite = G_malloc(sizeof(double));
*stats->overwrite =
- (double)G_get_raster_value_d(cell, IN_TYPE);
+ (double)Rast_get_d_value(cell, IN_TYPE);
}
else {
/* only add non-null cells to stats */
- cell_value = (double)G_get_raster_value_d(cell, IN_TYPE);
+ cell_value = (double)Rast_get_d_value(cell, IN_TYPE);
/* only add if within neighborhood */
if (WEIGHTS[i][j] != -1.0) {
/* get data needed for chosen statistic */
@@ -522,8 +522,7 @@
cell_input += CELL_IN_SIZE;
/* get next row from disk */
- G_get_raster_row(file_desc, cell_input, current_row + DATA_HEIGHT,
- IN_TYPE);
+ Rast_get_row(file_desc, cell_input, current_row + DATA_HEIGHT, IN_TYPE);
/* re-assign all row handlers below current replacement position */
j = PADDING_HEIGHT;
@@ -568,7 +567,7 @@
if (stats->num_values < min_cells) {
SET_NULL(cell_output, 1);
if (write_err)
- G_set_f_null_value(err_output, 1);
+ Rast_set_f_null_value(err_output, 1);
}
else {
if (stats->overwrite != NULL) {
@@ -583,10 +582,10 @@
}
/* write error/uncertainty output map? */
if (write_err) {
- G_set_raster_value_f(err_output,
- (FCELL) 1.0 -
- (stats->certainty / SUM_WEIGHTS),
- FCELL_TYPE);
+ Rast_set_f_value(err_output,
+ (FCELL) 1.0 -
+ (stats->certainty / SUM_WEIGHTS),
+ FCELL_TYPE);
}
}
/* advance cell pointers by one cell size */
@@ -1010,25 +1009,21 @@
}
/* open raster input map and get its storage type */
- mapset = G_find_cell2(input, "");
+ mapset = G_find_raster(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);
+ Rast_get_cellhd(input, mapset, &cellhd);
+ in_fd = Rast_open_old(input, mapset);
+ IN_TYPE = Rast_get_map_type(in_fd);
/* minimum and maximum values for interpolating range */
if (IN_TYPE == CELL_TYPE) {
- G_read_range(input, mapset, &int_range);
+ Rast_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);
+ Rast_read_fp_range(input, mapset, &fp_range);
min = (double)fp_range.min;
max = (double)fp_range.max;
}
@@ -1074,26 +1069,26 @@
}
}
}
- char data_type_string_in[21];
- char data_type_string_out[21];
+ char *data_type_string_in;
+ char *data_type_string_out;
if (IN_TYPE == CELL_TYPE) {
- G_strncpy(data_type_string_in, "integer", 20);
+ data_type_string_in = "integer";
}
else if (IN_TYPE == FCELL_TYPE) {
- G_strncpy(data_type_string_in, "single", 20);
+ data_type_string_in = "single";
}
else if (IN_TYPE == DCELL_TYPE) {
- G_strncpy(data_type_string_in, "double", 20);
+ data_type_string_in = "double";
}
if (OUT_TYPE == CELL_TYPE) {
- G_strncpy(data_type_string_out, "integer", 20);
+ data_type_string_out = "integer";
}
else if (OUT_TYPE == FCELL_TYPE) {
- G_strncpy(data_type_string_out, "single", 20);
+ data_type_string_out = "single";
}
else if (OUT_TYPE == DCELL_TYPE) {
- G_strncpy(data_type_string_out, "double", 20);
+ data_type_string_out = "double";
}
/* initialize data type dependent cell handling functions */
@@ -1123,7 +1118,8 @@
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);
+ Rast_set_null_value(CELL_INPUT[i], cols + (PADDING_WIDTH * 2),
+ IN_TYPE);
}
/*
@@ -1201,14 +1197,14 @@
*/
/* Open output map with right data type */
- out_fd = G_open_raster_new(output, OUT_TYPE);
+ out_fd = Rast_open_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);
+ CELL_OUTPUT = Rast_allocate_buf(OUT_TYPE);
/* initialize output row */
SET_NULL(CELL_OUTPUT, cols);
@@ -1216,14 +1212,14 @@
/* 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);
+ err_fd = Rast_open_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);
+ ERR_OUTPUT = Rast_allocate_buf(FCELL_TYPE);
/* initialize output row */
- G_set_f_null_value(ERR_OUTPUT, cols);
+ Rast_set_f_null_value(ERR_OUTPUT, cols);
}
/* row indices to handle input data buffer */
@@ -1250,16 +1246,16 @@
for (j = 0; j < PADDING_WIDTH; j++) {
cell_input += CELL_IN_SIZE;
}
- G_get_raster_row(in_fd, cell_input, i, IN_TYPE);
+ Rast_get_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);
+ Rast_put_row(out_fd, CELL_OUTPUT, OUT_TYPE);
if (parm.error->answer)
- G_put_raster_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
+ Rast_put_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
G_percent(current_row + 1, rows, 2);
current_row++;
}
@@ -1273,9 +1269,9 @@
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);
+ Rast_put_row(out_fd, CELL_OUTPUT, OUT_TYPE);
if (parm.error->answer)
- G_put_raster_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
+ Rast_put_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
G_percent(current_row + 1, rows, 2);
current_row++;
}
@@ -1285,25 +1281,25 @@
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);
+ Rast_get_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);
+ Rast_put_row(out_fd, CELL_OUTPUT, OUT_TYPE);
if (parm.error->answer)
- G_put_raster_row(err_fd, ERR_OUTPUT, FCELL_TYPE);
+ Rast_put_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);
+ Rast_close(out_fd);
+ Rast_close(in_fd);
if (parm.error->answer) {
- G_close_cell(err_fd);
+ Rast_close(err_fd);
}
/* Free memory */
@@ -1330,31 +1326,34 @@
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;
+ Rast_short_history(parm.output->answer, "raster", &hist);
+ Rast_put_cell_title(parm.output->answer,
+ "Result of interpolation/gap filling");
if (parm.dist_m->answer) {
- sprintf(hist.edhist[0],
- "Settings: mode=%s, distance (map units)=%.6f, power=%.3f",
- parm.mode->answer, radius, power);
+ Rast_append_format_history(&hist,
+ "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);
+ Rast_append_format_history(&hist,
+ "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);
+ Rast_append_format_history(&hist,
+ " min=%.3f, max=%.3f, min. points=%lu",
+ min, max, min_cells);
+ Rast_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);
+ Rast_short_history(parm.error->answer, "raster", &hist);
+ Rast_put_cell_title(parm.error->answer,
+ "Uncertainty of interpolation/gap filling");
+ Rast_append_format_history(&hist, "Result map: %s",
+ parm.output->answer);
+ Rast_append_format_history(&hist,
+ "Theoretic range is '0' (lowest) to '1' (highest).");
+ Rast_write_history(parm.error->answer, &hist);
}
finish = time(NULL);
Deleted: grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.01.png
===================================================================
(Binary files differ)
Deleted: grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.02.png
===================================================================
(Binary files differ)
Copied: grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.html (from rev 70879, grass-addons/grass7/raster/r.fill.gaps/description.html)
===================================================================
--- grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.html (rev 0)
+++ grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.html 2017-04-14 22:37:29 UTC (rev 70880)
@@ -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>.
Copied: grass-addons/grass7/raster/r.fill.gaps/r_fill_gaps_01.png (from rev 70879, grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.01.png)
===================================================================
(Binary files differ)
Copied: grass-addons/grass7/raster/r.fill.gaps/r_fill_gaps_02.png (from rev 70879, grass-addons/grass7/raster/r.fill.gaps/r.fill.gaps.02.png)
===================================================================
(Binary files differ)
More information about the grass-commit
mailing list