r.object.geometry: move from addons (#1918)

raster/Makefile

@@ -49,6 +49,7 @@ SUBDIRS = \
 	r.mode \
 	r.neighbors \
 	r.null \
+	r.object.geometry \
 	r.out.ascii \
 	r.out.bin \
 	r.out.gdal \

raster/r.object.geometry/Makefile

@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.object.geometry
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd

raster/r.object.geometry/main.c

@@ -0,0 +1,354 @@
+
+/****************************************************************************
+ *
+ * MODULE:       r.object.geometry
+ *
+ * AUTHOR(S):    Moritz Lennert
+ *               Markus Metz
+ *
+ * PURPOSE:      Fetch object geometry parameters.
+ *
+ * COPYRIGHT:    (C) 2016-2021 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ ***************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+/* compare two cell values
+ * return 0 if equal, 1 if different */
+static int cmp_cells(CELL a, CELL b, int a_null, int b_null)
+{
+    return (a_null + b_null == 1 || (a_null + b_null == 0 && a != b));
+}
+
+int main(int argc, char *argv[])
+{
+    int row, col, nrows, ncols;
+
+    struct Range range;
+    CELL min, max;
+    int in_fd;
+    int i;
+    struct GModule *module;
+    struct Option *opt_in;
+    struct Option *opt_out;
+    struct Option *opt_sep;
+    struct Flag *flag_m;
+    char *sep;
+    FILE *out_fp;
+    CELL *prev_in, *cur_in, *temp_in;
+    CELL cur, top, left;
+    int cur_null, top_null, left_null;
+    int len;
+    struct obj_geo
+    {
+        double area, perimeter, x, y;
+        int min_row, max_row, min_col, max_col; /* bounding box */
+        int num;
+    } *obj_geos;
+    double unit_area;
+    int n_objects;
+    int planimetric, compute_areas;
+    struct Cell_head cellhd;
+
+    G_gisinit(argv[0]);
+
+    /* Define the different options */
+
+    module = G_define_module();
+    G_add_keyword(_("raster"));
+    G_add_keyword(_("statistics"));
+    G_add_keyword(_("reclass"));
+    G_add_keyword(_("clumps"));
+    module->description =
+        _("Calculates geometry parameters for raster objects.");
+
+    opt_in = G_define_standard_option(G_OPT_R_INPUT);
+
+    opt_out = G_define_standard_option(G_OPT_F_OUTPUT);
+    opt_out->required = NO;
+
+    opt_sep = G_define_standard_option(G_OPT_F_SEP);
+
+    flag_m = G_define_flag();
+    flag_m->key = 'm';
+    flag_m->label = _("Use meters as units instead of cells");
+
+    /* parse options */
+    if (G_parser(argc, argv))
+        exit(EXIT_FAILURE);
+
+    sep = G_option_to_separator(opt_sep);
+    in_fd = Rast_open_old(opt_in->answer, "");
+
+    if (Rast_get_map_type(in_fd) != CELL_TYPE)
+        G_fatal_error(_("Input raster mus be of type CELL"));
+
+    if (opt_out->answer != NULL && strcmp(opt_out->answer, "-") != 0) {
+        if (!(out_fp = fopen(opt_out->answer, "w")))
+            G_fatal_error(_("Unable to open file <%s> for writing"),
+                          opt_out->answer);
+    }
+    else {
+        out_fp = stdout;
+    }
+
+    Rast_get_window(&cellhd);
+    nrows = cellhd.rows;
+    ncols = cellhd.cols;
+
+    /* allocate CELL buffers two columns larger than current window */
+    len = (ncols + 2) * sizeof(CELL);
+    prev_in = (CELL *) G_malloc(len);
+    cur_in = (CELL *) G_malloc(len);
+
+    /* fake a previous row which is all NULL */
+    Rast_set_c_null_value(prev_in, ncols + 2);
+
+    /* set left and right edge to NULL */
+    Rast_set_c_null_value(&cur_in[0], 1);
+    Rast_set_c_null_value(&cur_in[ncols + 1], 1);
+
+    Rast_read_range(opt_in->answer, "", &range);
+    Rast_get_range_min_max(&range, &min, &max);
+    if (Rast_is_c_null_value(&min) || Rast_is_c_null_value(&max))
+        G_fatal_error(_("Empty input map <%s>"), opt_in->answer);
+
+    /* REMARK: The following is only true if object ids are continuously numbered */
+    n_objects = max - min + 1;
+    obj_geos = G_malloc(n_objects * sizeof(struct obj_geo));
+
+    for (i = 0; i < n_objects; i++) {
+        obj_geos[i].area = 0;
+        obj_geos[i].perimeter = 0;
+        obj_geos[i].min_row = nrows;
+        obj_geos[i].max_row = 0;
+        obj_geos[i].min_col = ncols;
+        obj_geos[i].max_col = 0;
+        obj_geos[i].x = 0;
+        obj_geos[i].y = 0;
+        obj_geos[i].num = 0;
+    }
+
+    unit_area = 0.0;
+    if (flag_m->answer) {
+        switch (G_begin_cell_area_calculations()) {
+        case 0:                /* areas don't make sense, but ignore this for now */
+        case 1:
+            planimetric = 1;
+            unit_area = G_area_of_cell_at_row(0);
+            break;
+        default:
+            planimetric = 0;
+            break;
+        }
+    }
+    compute_areas = flag_m->answer && !planimetric;
+    G_begin_distance_calculations();
+
+    G_message(_("Calculating statistics"));
+    for (row = 0; row < nrows; row++) {
+        G_percent(row, nrows + 1, 2);
+
+        Rast_get_c_row(in_fd, cur_in + 1, row);
+        cur = cur_in[0];
+        cur_null = 1;
+        if (compute_areas)
+            unit_area = G_area_of_cell_at_row(row);
+        for (col = 1; col <= ncols; col++) {
+            left = cur;
+            cur = cur_in[col];
+            top = prev_in[col];
+
+            left_null = cur_null;
+            cur_null = Rast_is_c_null_value(&cur);
+            top_null = Rast_is_c_null_value(&top);
+
+            if (!cur_null) {
+                if (flag_m->answer) {
+                    obj_geos[cur - min].area += unit_area;
+                    obj_geos[cur - min].num += 1;
+                }
+                else {
+                    obj_geos[cur - min].area += 1;
+                }
+                obj_geos[cur - min].x +=
+                    Rast_col_to_easting(col - 0.5, &cellhd);
+                obj_geos[cur - min].y +=
+                    Rast_row_to_northing(row + 0.5, &cellhd);
+                if (obj_geos[cur - min].min_row > row)
+                    obj_geos[cur - min].min_row = row;
+                if (obj_geos[cur - min].max_row < row + 1)
+                    obj_geos[cur - min].max_row = row + 1;
+                if (obj_geos[cur - min].min_col > col)
+                    obj_geos[cur - min].min_col = col;
+                if (obj_geos[cur - min].max_col < col + 1)
+                    obj_geos[cur - min].max_col = col + 1;
+            }
+
+            if (cmp_cells(cur, top, cur_null, top_null)) {
+                if (flag_m->answer) {
+                    double perimeter;
+
+                    /* could be optimized */
+                    perimeter = G_distance(cellhd.west + col * cellhd.ew_res,
+                                           Rast_row_to_northing(row, &cellhd),
+                                           cellhd.west + (col +
+                                                          1) * cellhd.ew_res,
+                                           Rast_row_to_northing(row,
+                                                                &cellhd));
+
+                    if (!cur_null)
+                        obj_geos[cur - min].perimeter += perimeter;
+                    if (!top_null)
+                        obj_geos[top - min].perimeter += perimeter;
+                }
+                else {
+                    if (!cur_null)
+                        obj_geos[cur - min].perimeter += 1;
+                    if (!top_null)
+                        obj_geos[top - min].perimeter += 1;
+                }
+            }
+            if (cmp_cells(cur, left, cur_null, left_null)) {
+                if (flag_m->answer) {
+                    double perimeter;
+
+                    /* could be optimized */
+                    perimeter = G_distance(cellhd.west + col * cellhd.ew_res,
+                                           Rast_row_to_northing(row, &cellhd),
+                                           cellhd.west +
+                                           (col) * cellhd.ew_res,
+                                           Rast_row_to_northing(row + 1,
+                                                                &cellhd));
+
+                    if (!cur_null)
+                        obj_geos[cur - min].perimeter += perimeter;
+                    if (!left_null)
+                        obj_geos[left - min].perimeter += perimeter;
+                }
+                else {
+                    if (!cur_null)
+                        obj_geos[cur - min].perimeter += 1;
+                    if (!left_null)
+                        obj_geos[left - min].perimeter += 1;
+                }
+            }
+        }
+
+        /* last col, right borders */
+        if (flag_m->answer) {
+            double perimeter;
+
+            perimeter = G_distance(cellhd.east,
+                                   Rast_row_to_northing(row, &cellhd),
+                                   cellhd.east,
+                                   Rast_row_to_northing(row + 1, &cellhd));
+            if (!cur_null)
+                obj_geos[cur - min].perimeter += perimeter;
+        }
+        else {
+            if (!cur_null)
+                obj_geos[cur - min].perimeter += 1;
+        }
+
+        /* switch the buffers so that the current buffer becomes the previous */
+        temp_in = cur_in;
+        cur_in = prev_in;
+        prev_in = temp_in;
+    }
+
+    /* last row, bottom borders */
+    G_percent(row, nrows + 1, 2);
+    for (col = 1; col <= ncols; col++) {
+        top = prev_in[col];
+        top_null = Rast_is_c_null_value(&top);
+
+        if (flag_m->answer) {
+            double perimeter;
+
+            /* could be optimized */
+            perimeter = G_distance(cellhd.west + col * cellhd.ew_res,
+                                   Rast_row_to_northing(row, &cellhd),
+                                   cellhd.west + (col + 1) * cellhd.ew_res,
+                                   Rast_row_to_northing(row, &cellhd));
+
+            if (!top_null)
+                obj_geos[top - min].perimeter += perimeter;
+        }
+        else {
+            if (!top_null)
+                obj_geos[top - min].perimeter += 1;
+        }
+    }
+    G_percent(1, 1, 1);
+
+    Rast_close(in_fd);
+    G_free(cur_in);
+    G_free(prev_in);
+
+    G_message(_("Writing output"));
+    /* print table */
+    fprintf(out_fp, "cat%s", sep);
+    fprintf(out_fp, "area%s", sep);
+    fprintf(out_fp, "perimeter%s", sep);
+    fprintf(out_fp, "compact_square%s", sep);
+    fprintf(out_fp, "compact_circle%s", sep);
+    fprintf(out_fp, "fd%s", sep);
+    fprintf(out_fp, "mean_x%s", sep);
+    fprintf(out_fp, "mean_y");
+    fprintf(out_fp, "\n");
+
+    /* print table body */
+    for (i = 0; i < n_objects; i++) {
+        G_percent(i, n_objects - 1, 1);
+        /* skip empty objects */
+        if (obj_geos[i].area == 0)
+            continue;
+
+        fprintf(out_fp, "%d%s", min + i, sep);
+        fprintf(out_fp, "%f%s", obj_geos[i].area, sep);
+        fprintf(out_fp, "%f%s", obj_geos[i].perimeter, sep);
+        fprintf(out_fp, "%f%s",
+                4 * sqrt(obj_geos[i].area) / obj_geos[i].perimeter, sep);
+        fprintf(out_fp, "%f%s",
+                obj_geos[i].perimeter / (2 * sqrt(M_PI * obj_geos[i].area)),
+                sep);
+        /* log 1 = 0, so avoid that by always adding 0.001 to the area: */
+        fprintf(out_fp, "%f%s",
+                2 * log(obj_geos[i].perimeter) / log(obj_geos[i].area +
+                                                     0.001), sep);
+        if (!flag_m->answer)
+            obj_geos[i].num = obj_geos[i].area;
+        fprintf(out_fp, "%f%s", obj_geos[i].x / obj_geos[i].num, sep);
+        fprintf(out_fp, "%f", obj_geos[i].y / obj_geos[i].num);
+        /* object id: i + min */
+
+        /* TODO */
+        /* smoothness */
+        /* perimeter of bounding box / perimeter -> smoother objects have a higher smoothness value
+         * smoothness is in the range 0 < smoothness <= 1 */
+
+        /* bounding box perimeter */
+
+        /* bounding box size */
+
+        /* variance of X and Y to approximate bounding ellipsoid */
+
+        fprintf(out_fp, "\n");
+    }
+    if (out_fp != stdout)
+        fclose(out_fp);
+
+    exit(EXIT_SUCCESS);
+}

raster/r.object.geometry/r.object.geometry.html

@@ -0,0 +1,53 @@
+<h2>DESCRIPTION</h2>
+
+<p>
+<em>r.object.geometry</em> calculates form statistics of raster objects
+in the <b>input</b> map and writes it to the <b>output</b> text file
+ (or standard output if no output filename or '-' is given),
+with fields separated by the chosen <b>separator</b>.  Objects are defined
+ as clumps of adjacent cells with the same category value (e.g. output of 
+<em><a href="r.clump.html">r.clump</a></em> or
+<em><a href="i.segment.html">i.segment</a></em>).
+
+<p>
+By default, values are in pixels. If values in meters is desired, the user
+can set the <b>-m</b> flag. If the current working region is in lat-long or 
+has non-square pixels, using meters is recommended.
+
+<p>
+Statistics currently calculated are exactly the same as in
+<em><a href="v.to.db.html">v.to.db</a></em> (except for compact_square and 
+mean coordinates):
+
+<ul>
+<li>area</li>
+<li>perimeter</li>
+<li>compact_square (compactness compared to a square: 
+  <tt>compact_square = 4 * sqrt(area) / perimeter</tt>)
+<li>compact_circle (compactness compared to a circle: 
+  <tt>compact_circle = perimeter / ( 2 * sqrt(PI * area) )</tt>)</li>
+<li>fractal dimension ( <tt>fd = 2 * ( log(perimeter) / log(area + 0.001) )</tt> )</li>
+<li>mean x coordinate of object (in map units)</li>
+<li>mean y coordinate of object (in map units)</li>
+</ul>
+
+<h2>EXAMPLE</h2>
+
+<div class="code"><pre>
+g.region raster=soilsID
+r.object.geometry input=soilsID output=soils_geom.txt
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.segment.html">i.segment</a>,
+<a href="r.clump.html">r.clump</a>,
+<a href="v.to.db.html">v.to.db</a>
+</em>
+
+<h2>AUTHORS</h2>
+
+Moritz Lennert<br>
+Markus Metz (diagonal clump tracing)
+

raster/r.object.geometry/testsuite/data/file_meter.csv

@@ -0,0 +1,4 @@
+cat|area|perimeter|compact_square|compact_circle|fd|mean_x|mean_y
+1|750000000.000000|110000.000000|0.995859|1.133071|1.136081|625000.000000|237500.000000
+2|1500000000.000000|160000.000000|0.968246|1.165385|1.134278|655000.000000|225000.000000
+3|750000000.000000|110000.000000|0.995859|1.133071|1.136081|625000.000000|212500.000000

raster/r.object.geometry/testsuite/data/file_pixel.csv

@@ -0,0 +1,4 @@
+cat|area|perimeter|compact_square|compact_circle|fd|mean_x|mean_y
+1|4.000000|8.000000|1.000000|1.128379|2.999459|625000.000000|237500.000000
+2|8.000000|12.000000|0.942809|1.196827|2.389831|655000.000000|225000.000000
+3|4.000000|8.000000|1.000000|1.128379|2.999459|625000.000000|212500.000000

raster/r.object.geometry/testsuite/r_object_geometry_test.py

@@ -0,0 +1,122 @@
+"""
+Name:      r_object_geometry_test
+Purpose:   This script is to demonstrate a unit test for r.object.geometry
+           module.
+"""
+
+import os
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+testraster1 = """\
+north:   250000
+south:   200000
+east:    670000
+west:    610000
+rows:    4
+cols:    4
+
+1 1 2 2
+1 1 2 2
+3 3 2 2
+3 3 2 2
+"""
+
+
+class TestObjectGeometryPixel(TestCase):
+    """Test case for object geometry module"""
+
+    # Setup variables to be used for outputs
+    test_objects1 = "test_objects1"
+    output_file_pixel = "output_file_pixel.csv"
+
+    @classmethod
+    def setUpClass(cls):
+        """Imports test raster(s), ensures expected computational region and setup"""
+        cls.runModule(
+            "r.in.ascii", input="-", stdin_=testraster1, output=cls.test_objects1
+        )
+        cls.use_temp_region()
+        cls.runModule("g.region", raster=cls.test_objects1)
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region"""
+        cls.del_temp_region()
+
+    def tearDown(self):
+        """Remove the outputs created from the object geometry module
+
+        This is executed after each test run.
+        """
+        if os.path.isfile(self.output_file_pixel):
+            os.remove(self.output_file_pixel)
+        self.runModule("g.remove", flags="f", type="raster", name=self.test_objects1)
+
+    def test_object_geometry_pixel(self):
+        """Test to see if the outputs are created and are correct in pixel units"""
+        # run the object geometry module
+        self.assertModule(
+            "r.object.geometry", input=self.test_objects1, output=self.output_file_pixel
+        )
+        # check to see if output file exists
+        self.assertFileExists(self.output_file_pixel, msg="Output file does not exist")
+        # check if the output file is equal to the reference file
+        self.assertFilesEqualMd5(
+            self.output_file_pixel,
+            "data/file_pixel.csv",
+            msg="Output file is not equal to reference file",
+        )
+
+
+class TestObjectGeometryMeter(TestCase):
+    """Test case for object geometry module"""
+
+    # Setup variables to be used for outputs
+    test_objects1 = "test_objects1"
+    output_file_meter = "output_file_meter.csv"
+
+    @classmethod
+    def setUpClass(cls):
+        """Imports test raster(s), ensures expected computational region and setup"""
+        cls.runModule(
+            "r.in.ascii", input="-", stdin_=testraster1, output=cls.test_objects1
+        )
+        cls.use_temp_region()
+        cls.runModule("g.region", raster=cls.test_objects1)
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the temporary region"""
+        cls.del_temp_region()
+
+    def tearDown(self):
+        """Remove the outputs created from the object geometry module
+
+        This is executed after each test run.
+        """
+        if os.path.isfile(self.output_file_meter):
+            os.remove(self.output_file_meter)
+        self.runModule("g.remove", flags="f", type="raster", name=self.test_objects1)
+
+    def test_object_geometry_meter(self):
+        """Test to see if the outputs are created and are correct in meter units"""
+        # run the object geometry module
+        self.assertModule(
+            "r.object.geometry",
+            input=self.test_objects1,
+            output=self.output_file_meter,
+            flags="m",
+        )
+        # check to see if output file exists
+        self.assertFileExists(self.output_file_meter, msg="Output file does not exist")
+        # check if the output file is equal to the reference file
+        self.assertFilesEqualMd5(
+            self.output_file_meter,
+            "data/file_meter.csv",
+            msg="Output file is not equal to reference file",
+        )
+
+
+if __name__ == "__main__":
+    test()