-If the -c flag is used, r.reclass.area will skip the -creation of a clumped raster and assume that the input raster is +r.reclass.area removes areas (pixel clumps) from a raster +map that are smaller than a user specified lower area +size and/or greater than a user-specified upper +threshold (both in hectares). + +r.reclass.area provides two methods for +filtering: + +
-r.reclass.area input=zipcodes output=zipcodes_larger2000ha mode=greater value=2000 +r.reclass.area input=zipcodes output=zipcodes_larger2000ha greater=2000 r.report zipcodes_larger2000ha unit=h
r.reclass.area input=zipcodes output=zipcodes_minor1000ha mode=lesser value=1000 method=rmarea @@ -47,12 +66,16 @@EXAMPLES
SEE ALSO
-r.reclass, r.clump, -r.stats +r.reclass, +r.stats, +r.stats, +r.stats, +r.statsAUTHORS
NRCS,
-Markus Neteler +Markus Neteler,
+Stefan Blumentath diff --git a/scripts/r.reclass.area/r.reclass.area.md b/scripts/r.reclass.area/r.reclass.area.md index a20ff64d596..b46e7d4450c 100644 --- a/scripts/r.reclass.area/r.reclass.area.md +++ b/scripts/r.reclass.area/r.reclass.area.md @@ -1,10 +1,25 @@ ## DESCRIPTION -*r.reclass.area* reclasses a raster map greater or less than a user -specified area size (in hectares). +*r.reclass.area* removes areas (pixel clumps) from a raster map +that are smaller than a user specified **lower** area size and/or +greater than a user-specified **upper** threshold (both in hectares). -If the **-c** flag is used, *r.reclass.area* will skip the creation of a -clumped raster and assume that the input raster is already clumped. +*r.reclass.area* provides two **method**s for filtering: + +- *reclass* - raster-based filtering +- *rmarea* - vector-based filtering + +With the *reclass* method, the input raster map is clumped and +then reclassified using the clump size. With the **-m** flag, +*reclass* method will remove small areas already during clumping, +using the *minsize* option there. If the **-c** flag is used, +*r.reclass.area* will skip the creation of a clumped raster and +assume that the input raster map is already clumped. + +The *rmarea* method converts the raster map to vector and removes +areas outside the given thresholds using [v.clean](v.clean.md) and +[v.edit](v.edit.md). With the **-v** flag the **output** is the +filtered vector map instead of a raster. ## EXAMPLES @@ -19,7 +34,7 @@ r.report zipcodes unit=h Extract only areas greater than 2000 ha, NULL otherwise: ```sh -r.reclass.area input=zipcodes output=zipcodes_larger2000ha mode=greater value=2000 +r.reclass.area input=zipcodes output=zipcodes_larger2000ha greater=2000 r.report zipcodes_larger2000ha unit=h ``` @@ -30,10 +45,10 @@ r.report zipcodes_larger2000ha unit=h In this example, the ZIP code map in the North Carolina sample dataset is filtered for smaller areas which are substituted with the value of the respective adjacent area with largest shared boundary. Reclass by -substitutional removing of areas minor of 1000 ha: +substitutional removing of areas smaller than 1000 ha: ```sh -r.reclass.area input=zipcodes output=zipcodes_minor1000ha mode=lesser value=1000 method=rmarea +r.reclass.area input=zipcodes output=zipcodes_minor1000ha lesser=1000 method=rmarea ```  @@ -41,10 +56,12 @@ r.reclass.area input=zipcodes output=zipcodes_minor1000ha mode=lesser value=1000 ## SEE ALSO -*[r.reclass](r.reclass.md), [r.clump](r.clump.md), -[r.stats](r.stats.md)* +*[r.clump](r.clump.md), [r.reclass](r.reclass.md), +[r.stats](r.stats.md), [v.to.rast](v.to.rast.md), +[v.clean](v.clean.md), [v.edit](v.edit.md)* ## AUTHORS NRCS, Markus Neteler +Stefan Blumentath diff --git a/scripts/r.reclass.area/r.reclass.area.py b/scripts/r.reclass.area/r.reclass.area.py index ff43f26e454..c0546ef163a 100755 --- a/scripts/r.reclass.area/r.reclass.area.py +++ b/scripts/r.reclass.area/r.reclass.area.py @@ -7,13 +7,14 @@ # Converted to Python by Glynn Clements # Added rmarea method by Luca Delucchi # PURPOSE: Reclasses a raster map greater or less than user specified area size (in hectares) -# COPYRIGHT: (C) 1999,2008,2014 by the GRASS Development Team +# COPYRIGHT: (C) 1999-2026 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. # ############################################################################# +# 04/2026: rewrite with Tools API and new options (Stefan Blumentrath) # 10/2013: added option to use a pre-clumped input map (Eric Goddard) # 8/2012: added fp maps support, cleanup, removed tabs AK # 3/2007: added label support MN @@ -39,24 +40,40 @@ # %option # % key: value # % type: double -# % description: Value option that sets the area size limit (in hectares) -# % required: yes +# % description: Value option that sets the area size limit (in hectares) (deprecated) +# % required: no # % guisection: Area # %end # %option # % key: mode # % type: string -# % description: Lesser or greater than specified value +# % description: Keep only areas lesser or greater than specified value (deprecated) # % options: lesser,greater -# % required: yes +# % required: no +# % guisection: Area +# %end + +# %option +# % key: lower +# % type: double +# % description: Value for the lower area size limit (in hectares). Pixel clumps with areas below this value are removed. +# % required: no +# % guisection: Area +# %end + +# %option +# % key: upper +# % type: double +# % description: Value for the upper area size limit (in hectares). Pixel clumps with areas above this value are removed. +# % required: no # % guisection: Area # %end # %option # % key: method # % type: string -# % description: Method used for reclassification +# % description: Method used for area size filtering # % options: reclass,rmarea # % answer: reclass # % guisection: Area @@ -69,174 +86,300 @@ # %flag # % key: d -# % description: Clumps including diagonal neighbors +# % description: Generate clumps including diagonal neighbors # %end -import sys -import os +# %flag +# % key: i +# % description: Invert filter (remove clumps larger than lower and smaller than upper) +# %end + +# %flag +# % key: m +# % description: Apply lower area threshold as "minsize" during clumping +# %end + +# %flag +# % key: v +# % description: Output vector map, do not convert back to raster (only for method 'rmarea') +# %end + +# %rules +# % required: lower, upper, mode +# % requires: mode, value +# % requires: -m, lower +# % excludes: -c, -d, -m +# %end + + import atexit -import grass.script as gs -from grass.script.utils import decode, encode +import io +import sys -TMPRAST = [] +from functools import partial +from math import ceil, inf, prod +from typing import Literal +import grass.script as gs +from grass.tools import Tools -def reclass(inf, outf, lim, clump, diag, les): - infile = inf - outfile = outf - lesser = les - limit = lim - clumped = clump - diagonal = diag +TEMP_PREFIX = gs.tempname(12) - s = gs.read_command("g.region", flags="p") - s = decode(s) - kv = gs.parse_key_val(s, sep=":") - s = kv["projection"].strip().split() - if s == "0": - gs.fatal(_("xy-locations are not supported")) - gs.fatal(_("Need projected data with grids in meters")) - if not gs.find_file(infile)["name"]: - gs.fatal(_("Raster map <%s> not found") % infile) +def cleanup() -> None: + """Delete all temporary maps.""" + tools = Tools(capture_output=False) + tools.g_remove( + flags="fb", + type="all", + pattern=f"{TEMP_PREFIX}_*", + superquiet=True, + ) - if clumped and diagonal: - gs.fatal(_("flags c and d are mutually exclusive")) - if clumped: - clumpfile = infile - else: - clumpfile = "%s.clump.%s" % (infile.split("@")[0], outfile) - TMPRAST.append(clumpfile) - - if not gs.overwrite(): - if gs.find_file(clumpfile)["name"]: - gs.fatal(_("Temporary raster map <%s> exists") % clumpfile) - if diagonal: - gs.message( - _("Generating a clumped raster file including diagonal neighbors...") - ) - gs.run_command("r.clump", flags="d", input=infile, output=clumpfile) - else: - gs.message(_("Generating a clumped raster file ...")) - gs.run_command("r.clump", input=infile, output=clumpfile) - - if lesser: - gs.message( +def check_projection() -> tuple[float, float]: + """Check if current projection is supported. + + Issues a fatal message if the tool does not reliably + work in the current type of projection (XY, latlon) + """ + tools = Tools(capture_output=True) + proj = tools.g_region(flags="p", format="json").json + if proj["crs"]["type_code"] == "0": + gs.fatal( _( - "Generating a reclass map with area size less than " - "or equal to %f hectares..." - ) - % limit + "xy-locations are not supported\n" + "Need projected data with grids in meters", + ), ) + in_proj = gs.parse_key_val(tools.g_proj(flags="p", format="shell").text) + if in_proj["unit"].lower() == "degree": + gs.fatal(_("Latitude-longitude locations are not supported")) + if in_proj["name"].lower() == "xy_location_unprojected": + gs.fatal(_("xy-locations are not supported")) + return proj["ewres"], proj["nsres"] + + +def get_clumpfile( + input_map: str, + resolution: tuple[float, float], + minsize: float | None = None, + *, + clump_flags: Literal["d"] | None = None, +) -> str: + """Generate a clump-raster map.""" + tools = Tools(capture_output=False) + clump_map = f"{TEMP_PREFIX}_clump" + + if minsize: + minsize = ceil((minsize * 10000.0) / prod(resolution)) + + msg = _("Generating a clumped raster file") + if clump_flags: + msg += _(" including diagonal neighbors") + if minsize: + msg += _(" with minimum %s pixels") % minsize + gs.verbose(msg) + tools.r_clump( + flags=clump_flags, + input=input_map, + minsize=minsize, + output=clump_map, + quiet=True, + ) + return clump_map + + +def reclass( + clump_map: str, + output_map: str | None, + lower: float | None = None, + upper: float | None = None, + input_map: str | None = None, +) -> None: + """Perform raster-based filtering.""" + tools = Tools(capture_output=True) + stats_input = clump_map + expected_fields_number = 4 + verbose_message = _("Generating a reclass map with area size") + range_filter_message = "" + + # Define lower threshold + if lower: + range_filter_message += _(" larger than or equal to %f hectares...") % lower else: - gs.message( - _( - "Generating a reclass map with area size greater " - "than or equal to %f hectares..." - ) - % limit - ) + lower = -inf - recfile = outfile + ".recl" - TMPRAST.append(recfile) + # Define upper threshold + if upper: + if range_filter_message: + range_filter_message += " and" + range_filter_message += _(" smaller than or equal to %f hectares...") % upper + else: + upper = inf + gs.verbose(verbose_message + range_filter_message) + recfile = f"{TEMP_PREFIX}_recl" sflags = "aln" - if gs.raster_info(infile)["datatype"] in {"FCELL", "DCELL"}: - sflags += "i" - p1 = gs.pipe_command("r.stats", flags=sflags, input=(clumpfile, infile), sep=";") - p2 = gs.feed_command("r.reclass", input=clumpfile, output=recfile, rules="-") + if input_map: + # For non-null clumps r.stats should produce 5 columns + expected_fields_number = 5 + stats_input = (clump_map, input_map) + if gs.raster_info(input_map)["datatype"] in {"FCELL", "DCELL"}: + sflags += "i" + clump_areas = ( + tools.r_stats( + flags=sflags, + input=stats_input, + sep=";", + quiet=True, + ) + .text.rstrip() + .split("\n") + ) + rules = "" - for line in p1.stdout: - f = decode(line).rstrip(os.linesep).split(";") - if len(f) < 5: + for line in clump_areas: + f = line.rstrip().split(";") + if len(f) < expected_fields_number: continue - hectares = float(f[4]) * 0.0001 - test = hectares <= limit if lesser else hectares >= limit - if test: - rules += "%s = %s %s\n" % (f[0], f[2], f[3]) + hectares = float(f[-1]) * 0.0001 + if lower <= hectares <= upper: + rules += "{} = {} {}\n".format(*[f[i] for i in [0, 2, -2]]) if rules: - p2.stdin.write(encode(rules)) - p1.wait() - p2.stdin.close() - p2.wait() - if p2.returncode != 0: - if lesser: - gs.fatal( - _("No areas of size less than or equal to %f hectares found.") % limit - ) - else: - gs.fatal( - _("No areas of size greater than or equal to %f hectares found.") - % limit - ) - gs.mapcalc("$outfile = $recfile", outfile=outfile, recfile=recfile) - - -def rmarea(infile, outfile, thresh, coef): - # transform user input from hectares to map units (kept this for future) - # thresh = thresh * 10000.0 / (float(coef)**2) - # grass.debug("Threshold: %d, coeff linear: %s, coef squared: %d" % - # (thresh, coef, (float(coef)**2)), 0) - - # transform user input from hectares to meters because currently v.clean - # rmarea accept only meters as threshold - thresh *= 10000.0 - vectfile = "%s_vect_%s" % (infile.split("@")[0], outfile) - TMPRAST.append(vectfile) - gs.run_command("r.to.vect", input=infile, output=vectfile, type="area") - cleanfile = "%s_clean_%s" % (infile.split("@")[0], outfile) - TMPRAST.append(cleanfile) - gs.run_command( - "v.clean", input=vectfile, output=cleanfile, tool="rmarea", threshold=thresh - ) + tools.r_reclass( + input=clump_map, + output=recfile, + rules=io.StringIO(rules), + quiet=True, + ) + else: + fatal_message = _("No areas of size ") + gs.fatal(fatal_message + range_filter_message) + tools.r_mapcalc(expression=f"{output_map} = {recfile}") + + +def rmarea( + input_map: str, + output_map: str, + lower: float | None = None, + upper: float | None = None, + *, + return_vector: bool = False, +) -> None: + """Perform vector-based filtering.""" + tools = Tools(capture_output=False) + + if gs.raster_info(input_map)["datatype"] != "CELL": + gs.warning( + _("Raster map <%s> is not of type 'CELL'. Results may not be as expected."), + ) - gs.run_command( - "v.to.rast", input=cleanfile, output=outfile, use="attr", attrcolumn="value" + # Setup temporary map names with pattern + vectfile = f"{TEMP_PREFIX}_vect" + cleanfile = f"{TEMP_PREFIX}_clean" + edit_vector = vectfile + + # Convert input map to vector raster values as category + tools.r_to_vect( + flags="v", + input=input_map, + output=vectfile, + type="area", + quiet=True, ) + # Apply lower threshold + # transform user input from hectares to meters because currently v.clean + # rmarea accepts only meters as threshold + if lower: + lower *= 10000.0 + tools.v_clean( + input=vectfile, + output=cleanfile, + tool="rmarea", + threshold=lower, + quiet=True, + ) + edit_vector = cleanfile + + # Apply upper threshold + if upper: + tools.v_to_db( + map=edit_vector, + type="centroid", + option="area", + columns="area_ha", + units="hectares", + quiet=True, + ) + tools.v_edit( + map=edit_vector, + tool="delete", + where=f"area_ha >= {upper}", + superquiet=True, + ) -def main(): - infile = options["input"] - value = options["value"] - mode = options["mode"] - outfile = options["output"] - global method - method = options["method"] - clumped = flags["c"] - diagonal = flags["d"] - - # check for unsupported locations - in_proj = gs.parse_command("g.proj", flags="p", format="shell") - if in_proj["unit"].lower() == "degree": - gs.fatal(_("Latitude-longitude locations are not supported")) - if in_proj["name"].lower() == "xy_location_unprojected": - gs.fatal(_("xy-locations are not supported")) - - # check lesser and greater parameters - limit = float(value) - if mode == "greater" and method == "rmarea": - gs.fatal(_("You have to specify mode='lesser' with method='rmarea'")) + # Produce final output map (raster / vector) + if not return_vector: + tools.v_to_rast( + input=edit_vector, + output=output_map, + use="cat", + quiet=True, + ) + else: + tools.v_clean( + flags="c", + input=edit_vector, + output=output_map, + tool=["rmdangle", "rmbridge"], + threshold=-1, + quiet=True, + ) - if not gs.find_file(infile)["name"]: - gs.fatal(_("Raster map <%s> not found") % infile) - if method == "reclass": - reclass(infile, outfile, limit, clumped, diagonal, mode == "lesser") - elif method == "rmarea": - rmarea(infile, outfile, limit, in_proj["meters"]) +def main() -> None: + """Run area filtering.""" + if not gs.find_file(options["input"])["name"]: + gs.fatal(_("Raster map <%s> not found") % options["input"]) - gs.message(_("Generating output raster map <%s>...") % outfile) + # Check for unsupported locations + resolution = check_projection() + # Give deprecation warnings + if options["mode"]: + gs.warning( + _( + "Optionsand are deprecated. " + "Use 'lower' and/or 'upper'.", + ), + ) + options["lower" if options["mode"] == "greater" else "upper"] = options["value"] + + # Define filter thresholds + lower = float(options["lower"]) if options["lower"] else None + upper = float(options["upper"]) if options["upper"] else None + + # Compute clump map if needed + clump_map = ( + options["input"] + if flags["c"] or options["method"] == "rmarea" + else get_clumpfile( + options["input"], + resolution=resolution, + minsize=lower if flags["m"] else None, + clump_flags="d" if flags["d"] else None, + ) + ) -def cleanup(): - """!Delete temporary maps""" - TMPRAST.reverse() # reclassed map first - for mapp in TMPRAST: - if method == "rmarea": - gs.run_command("g.remove", flags="f", type="vector", name=mapp, quiet=True) - else: - gs.run_command("g.remove", flags="f", type="raster", name=mapp, quiet=True) + # Filter clumped map + filter_method = ( + partial(reclass, input_map=options["input"]) + if options["method"] == "reclass" + else partial(rmarea, return_vector=flags["v"]) + ) + filter_method(clump_map, options["output"], lower=lower, upper=upper) if __name__ == "__main__": diff --git a/scripts/r.reclass.area/testsuite/test_r_reclass_area.py b/scripts/r.reclass.area/testsuite/test_r_reclass_area.py index 96c922bda05..b2a59fdaee6 100644 --- a/scripts/r.reclass.area/testsuite/test_r_reclass_area.py +++ b/scripts/r.reclass.area/testsuite/test_r_reclass_area.py @@ -17,6 +17,9 @@ class TestReclassArea(TestCase): input = "geology_30m" output = "reclassarea" value = "20" + upper = "4000" + lower = "20" + map_list = [] @classmethod def setUpClass(cls): @@ -28,43 +31,198 @@ def tearDownClass(cls): cls.del_temp_region() cls.runModule( "g.remove", - type="raster", + type="all", flags="f", - name=(cls.output + "Greater", cls.output + "Lesser"), + name=cls.map_list, ) - def test_reclassaeaGreater(self): - """Testing r.reclass.area with greater""" + def test_reclass_greater(self) -> None: + """Testing r.reclass.area with method reclass and greater (old mode).""" + output_map = f"{self.output}_greater" + ref_max = 946 + ref_min = 217 self.assertModule( "r.reclass.area", input=self.input, - output=self.output + "Greater", + output=output_map, value=self.value, mode="greater", method="reclass", ) + self.map_list.append(output_map) self.assertRasterMinMax( - map=self.output + "Greater", - refmin=200, - refmax=1000, - msg="Range of data: min = 200 max = 1000", + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", ) - def test_reclassareaLesser(self): - """Testing r.reclass.area with lesser""" + def test_reclass_lower(self) -> None: + """Testing r.reclass.area with method reclass and lower.""" + output_map = f"{self.output}_lower" + ref_max = 1000 + ref_min = 200 self.assertModule( "r.reclass.area", input=self.input, - output=self.output + "Lesser", + output=output_map, + lower=self.lower, + method="reclass", + ) + self.map_list.append(output_map) + self.assertRasterMinMax( + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", + ) + + def test_reclass_lesser(self) -> None: + """Testing r.reclass.area with method reclass and lesser (old mode).""" + output_map = f"{self.output}_lesser" + ref_max = 1000 + ref_min = 900 + self.assertModule( + "r.reclass.area", + input=self.input, + output=output_map, value=self.value, mode="lesser", method="reclass", ) + self.map_list.append(output_map) + self.assertRasterMinMax( + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", + ) + + def test_reclass_upper(self) -> None: + """Testing r.reclass.area with method reclass and upper.""" + output_map = f"{self.output}_upper" + ref_max = 1000 + ref_min = 262 + self.assertModule( + "r.reclass.area", + input=self.input, + output=output_map, + upper=self.upper, + method="reclass", + ) + self.map_list.append(output_map) self.assertRasterMinMax( - map=self.output + "Lesser", - refmin=900, - refmax=1000, - msg="Range of data: min = 900 max = 1000", + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", + ) + + def test_reclass_lower_upper(self) -> None: + """Testing r.reclass.area with method reclass and lower and upper.""" + output_map = f"{self.output}_lower_upper" + ref_max = 1000 + ref_min = 200 + self.assertModule( + "r.reclass.area", + input=self.input, + output=output_map, + lower=self.lower, + upper=self.upper, + method="reclass", + ) + self.map_list.append(output_map) + self.assertRasterMinMax( + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", + ) + + def test_rmarea_lower(self) -> None: + """Testing r.reclass.area with method rmarea and lower.""" + output_map = f"{self.output}_rmarea_lower" + ref_max = 946 + ref_min = 217 + self.assertModule( + "r.reclass.area", + input=self.input, + output=output_map, + lower=self.lower, + method="rmarea", + ) + self.map_list.append(output_map) + self.assertRasterMinMax( + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", + ) + + def test_rmarea_upper(self) -> None: + """Testing r.reclass.area with method rmarea and upper.""" + output_map = f"{self.output}_rmarea_upper" + ref_max = 948 + ref_min = 262 + self.assertModule( + "r.reclass.area", + input=self.input, + output=output_map, + upper=self.upper, + method="rmarea", + ) + self.map_list.append(output_map) + self.assertRasterMinMax( + map=output_map, + refmin=ref_min, + refmax=ref_max, + msg=f"Range of data: min = {ref_min} max = {ref_max}", + ) + + def test_rmarea_lower_upper_vector(self) -> None: + """Testing r.reclass.area with method rmarea, lower and upper.""" + output_map = f"{self.output}_rmarea_lower_upper" + self.assertModule( + "r.reclass.area", + flags="vd", + input=self.input, + output=output_map, + lower=self.lower, + upper=self.upper, + method="rmarea", + ) + self.map_list.append(output_map) + self.assertVectorExists( + output_map, + msg=f"Output map {output_map} not found.", + ) + self.assertVectorFitsTopoInfo( + vector=output_map, + reference={"areas": 11, "primitives": 43}, + msg=f"Toplogy of output map {output_map} does not match reference.", + ) + + def test_rmarea_lower_upper_vector_basins(self) -> None: + """Testing r.reclass.area with method rmarea, lower and upper and basins input.""" + output_map = f"{self.output}_rmarea_lower_upper_basin" + self.assertModule( + "r.reclass.area", + flags="vd", + input="basin", + output=output_map, + lower=self.lower, + upper=self.upper, + method="rmarea", + ) + self.map_list.append(output_map) + self.assertVectorExists( + output_map, + msg=f"Output map {output_map} not found.", + ) + self.assertVectorFitsTopoInfo( + vector=output_map, + reference={"areas": 14, "primitives": 82}, + msg=f"Toplogy of output map {output_map} does not match reference.", ) diff --git a/scripts/r.reclass.area/testsuite/testrra.py b/scripts/r.reclass.area/testsuite/testrra.py index 5e11bf6aecc..c27e07f11b7 100644 --- a/scripts/r.reclass.area/testsuite/testrra.py +++ b/scripts/r.reclass.area/testsuite/testrra.py @@ -89,9 +89,9 @@ def test_method_rmarea(self): ) self.assertRasterMinMax( map=self.output, - refmin=27603, - refmax=27607, - msg="Output Map in degrees must be between 27603 and 27607", + refmin=27511, + refmax=27610, + msg="Output Map in degrees must be between 27511 and 27610", )