v.surf.icw: Python3 and modern GRASS compatibility, fix bug introduced in trac # 2574 #1200
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,8 +1,7 @@ | ||
| #!/usr/bin/env python3 | ||
| ############################################################################# | ||
| # | ||
| # MODULE: v.surf.icw | ||
| # | ||
| # AUTHOR: M. Hamish Bowman, Dunedin, New Zealand | ||
| # Originally written aboard the NZ DoC ship M/V Renown, | ||
|
|
@@ -13,7 +12,7 @@ | |
| # PURPOSE: Like IDW interpolation, but distance is cost to get to any | ||
| # other site. | ||
| # | ||
| # COPYRIGHT: (c) 2003-2024 Hamish Bowman | ||
| # This program is free software under the GNU General Public | ||
| # License (>=v2). Read the file COPYING that comes with GRASS | ||
| # for details. | ||
|
|
@@ -44,25 +43,14 @@ | |
| # % keyword: interpolation | ||
| # % keyword: ICW | ||
| # %End | ||
| # %option G_OPT_V_INPUT | ||
| # % label: Name of existing vector points map containing seed data | ||
| # %end | ||
| # %option G_OPT_DB_COLUMN | ||
| # % description: Column name in points map that contains data values | ||
| # % required : yes | ||
| # %end | ||
| # %option G_OPT_R_OUTPUT | ||
| # %end | ||
| # %option | ||
| # % key: cost_map | ||
|
|
@@ -79,19 +67,11 @@ | |
| # % options: 1-6 | ||
| # % required : no | ||
| # %end | ||
| # %option G_OPT_V_FIELD | ||
| # % type: integer | ||
| # % label: Layer number of data in points map | ||
| # %end | ||
| # %option G_OPT_DB_WHERE | ||
| # %end | ||
|
|
||
| ##%option | ||
|
|
@@ -112,11 +92,8 @@ | |
| # % key: r | ||
| # % description: Use (d^n)*log(d) instead of 1/(d^n) for radial basis function | ||
| # %end | ||
| # %option G_OPT_M_NPROCS | ||
| # % options: 1-1024 | ||
| # % description: Number of parallel processes to launch | ||
| # %end | ||
|
|
||
|
|
@@ -135,14 +112,15 @@ | |
| def cleanup(): | ||
| grass.verbose(_("Cleanup..")) | ||
| tmp_base = "tmp_icw_" + str(os.getpid()) + "_" | ||
| # grass.run_command('g.list', type='raster', mapset='.', flags='p') | ||
| result = grass.list_strings("raster", pattern=tmp_base + "*", mapset=".") | ||
| if len(result) > 0: | ||
| grass.run_command( | ||
| "g.remove", flags="f", type="raster", pattern=tmp_base + "*", quiet=True | ||
| ) | ||
|
Comment on lines
+115
to
+120
There was a problem hiding this comment. Why the direct call to g.remove does not work without checking the existence first? Too verbose? There is |
||
|
|
||
|
|
||
| def main(): | ||
| pts_input = options["input"] | ||
| output = options["output"] | ||
| cost_map = options["cost_map"] | ||
|
|
@@ -151,7 +129,7 @@ def main(): | |
| friction = float(options["friction"]) | ||
| layer = options["layer"] | ||
| where = options["where"] | ||
| workers = int(options["nprocs"]) | ||
|
|
||
| if workers == 1 and "WORKERS" in os.environ: | ||
| workers = int(os.environ["WORKERS"]) | ||
|
|
@@ -267,20 +245,21 @@ def main(): | |
| northing = position[1] | ||
| cat = int(position[-1]) | ||
|
|
||
| # FIXME: layer=layer probably needed here | ||
| # retrieve data value from vector's attribute table: | ||
| data_value = grass.vector_db_select(pts_input, columns=column)["values"][cat][0] | ||
|
|
||
| if not data_value: | ||
| grass.verbose( | ||
| _("Site %d of %d, e=%.4f n=%.4f cat=%d data=?") | ||
| % (num, n, float(easting), float(northing), cat) | ||
| ) | ||
| grass.verbose(_(" -- Skipping, no data here.")) | ||
| del points_list[num - 1] | ||
| n -= 1 | ||
| continue | ||
| else: | ||
| grass.verbose( | ||
| _("Site %d of %d, e=%.4f n=%.4f cat=%d data=%.8g") | ||
| % (num, n, float(easting), float(northing), cat, float(data_value)) | ||
| ) | ||
|
|
@@ -296,7 +275,7 @@ def main(): | |
| .split("|")[-1] | ||
| ) | ||
| if rast_val == "*": | ||
| grass.verbose(_(" -- Skipping, point lays outside of cost_map.")) | ||
| del points_list[num - 1] | ||
| n -= 1 | ||
| continue | ||
|
|
@@ -327,6 +306,7 @@ def main(): | |
| if proc[i].wait() != 0: | ||
| grass.fatal(_("Problem running %s") % "r.cost") | ||
|
|
||
| grass.verbose("\n") | ||
| grass.message(_("Removing anomalies at site positions ...")) | ||
|
|
||
| proc = {} | ||
|
|
@@ -345,7 +325,7 @@ def main(): | |
| ) | ||
| # stall to wait for the nth worker to complete, | ||
| if (i + 1) % workers == 0: | ||
| # print('stalling ...') | ||
| proc[i].wait() | ||
|
|
||
| # make sure everyone is finished | ||
|
|
@@ -393,7 +373,7 @@ def main(): | |
| ) | ||
| # stall to wait for the nth worker to complete, | ||
| if (i + 1) % workers == 0: | ||
| # print('stalling ...') | ||
| proc[i].wait() | ||
|
|
||
| # r.patch in=1by_cost_site_sqrd.${NUM},tmp_idw_cost_val_$$ out=1by_cost_site_sqrd.${NUM} --o | ||
|
|
@@ -415,8 +395,7 @@ def main(): | |
|
|
||
| ####################################################### | ||
| #### Step 3) find sum(cost^2) | ||
| grass.verbose("\n" + _("Finding sum of squares ...")) | ||
|
|
||
| # todo: test if MASK exists already, fatal exit if it does? | ||
| if post_mask: | ||
|
|
@@ -425,14 +404,12 @@ def main(): | |
|
|
||
| grass.message(_("Summation of cost weights ...")) | ||
|
|
||
| TMP_FILE = grass.tempfile() | ||
| with open(TMP_FILE, "w") as maplist: | ||
| for i in range(1, n + 1): | ||
| mapname = "%s1by_cost_site_sq.%05d" % (tmp_base, i) | ||
| maplist.write(mapname + "\n") | ||
| maplist.close() | ||
|
|
||
| # grass.run_command('g.list', type = 'raster', mapset = '.') | ||
|
|
||
|
|
@@ -444,13 +421,15 @@ def main(): | |
| except CalledModuleError: | ||
| grass.fatal(_("Problem running %s") % "r.series") | ||
|
|
||
| grass.try_remove(TMP_FILE) | ||
|
|
||
| if post_mask: | ||
| grass.message(_("Removing post_mask <%s>"), post_mask) | ||
| grass.run_command("g.remove", flags="f", name="MASK", quiet=True) | ||
|
|
||
| ####################################################### | ||
| #### Step 4) ( 1/di^2 / sum(1/d^2) ) * ai | ||
| grass.verbose("\n") | ||
| grass.message(_("Creating partial weights ...")) | ||
|
|
||
| proc = {} | ||
|
|
@@ -459,17 +438,18 @@ def main(): | |
| easting = position[0] | ||
| northing = position[1] | ||
| cat = int(position[-1]) | ||
| # FIXME: layer=layer probably needed here | ||
| data_value = grass.vector_db_select(pts_input, columns=column)["values"][cat][0] | ||
| data_value = float(data_value) | ||
|
|
||
| # failsafe: at this point the data values should all be valid | ||
| if not data_value: | ||
| grass.verbose(_("Site %d of %d, cat = %d, data value = ?") % (num, n, cat)) | ||
| grass.verbose(_(" -- Skipping, no data here. [Probably programmer error]")) | ||
| n -= 1 | ||
| continue | ||
| else: | ||
| grass.verbose( | ||
| _("Site %d of %d, cat = %d, data value = %.8g") | ||
| % (num, n, cat, data_value) | ||
| ) | ||
|
|
@@ -485,7 +465,7 @@ def main(): | |
| .split("|")[-1] | ||
| ) | ||
| if rast_val == "*": | ||
| grass.verbose( | ||
| _( | ||
| " -- Skipping, point lays outside of cost_map. [Probably programmer error]" | ||
| ) | ||
|
|
@@ -534,30 +514,39 @@ def main(): | |
| # grass.run_command('g.list', type = 'raster', mapset = '.') | ||
|
|
||
| ####################################################### | ||
| grass.verbose("\n") | ||
| grass.message(_("Calculating final values ...")) | ||
|
|
||
| TMP_FILE = grass.tempfile() | ||
|
There was a problem hiding this comment. See https://github.com/OSGeo/grass/blob/main/doc/development/style_guide.md#temporary-files g.tempfile is meant for large data, in this case standard system tempfile seems better. There was a problem hiding this comment. See also how to deal with temporary maps in the guide. Right now you call cleanup function explicitly, but does it get called again since it's registered with atexit? There was a problem hiding this comment.
It gets called on exit a second time, but only tries to remove if something actually exists. It could be improved, tmp_base should be a global variable for one thing. There was a problem hiding this comment.
Any particular reason why? I don't really mind either way but g.tempfile only takes 1.36 milliseconds to complete for me, is only run twice during a long running module, doesn't pollute the filesystem outside of the mapset's .tmp dir, and gets a second chance to be cleared when the GRASS session ends/launches if the original core.try_remove() fails for some reason. g.tempfile just tries to create the file and reports back what its name is, it doesn't know or care about what ends up in it. For sure large files should use g.tempfile so that the disk space is used in MAPSET and we don't risk filling up a small /tmp partition. There was a problem hiding this comment. Yes, g.tempfile is meant for large files and also for Bash/shells, but otherwise there is little reason to use it. 1) In general, relying on a standard library, in this case Python standard library, is better than relying on our custom mechanism. The library provides standardized, well-document and well-tested tools with known limits. 2) Security concerns are known and addressed. 3) Code analyzers can understand the code better. 4) Other readers need to know only standard Python, not GRASS GIS specific API. 5) Creating files in the mapset can be significantly slower than creating them in system's tmp dir. That dir is meant to be fast, that's at least what most application would expect. /tmp is sometimes mounted to memory. On the other hand, mapset dir may be on a slow network drive for storage of GRASS projects, large data, or because of parallel, multi-node processing in the same GRASS project. Thinking about this more, some Python functions for tmp files and dirs allow you to pick the base tmp dir, so you can use the standard Python functions with mapset's .tmp if large files are a concern. (I used a custom dir approach to implement a temporary mapset session (TemporaryMapsetSession, create_temporary_mapset).) If that seems like a good way, we may want to document that g.tempfile is meant only for Bash/shells scripts and perhaps add a convenient wrapper for Python functions to work over mapset's tmp dir. Another take on the situation is that the standard Python functions (and related code quality checks) are designed in a way that the code should be sure everything is properly deleted at the right time. The |
||
| with open(TMP_FILE, "w") as maplist: | ||
| for i in range(1, n + 1): | ||
| mapname = "%spartial.%05d" % (tmp_base, i) | ||
| maplist.write(mapname + "\n") | ||
| maplist.close() | ||
|
|
||
| try: | ||
| grass.run_command("r.series", method="sum", file=TMP_FILE, output=output) | ||
| except CalledModuleError: | ||
| grass.fatal(_("Problem running %s") % "r.series") | ||
|
|
||
| grass.try_remove(TMP_FILE) | ||
|
|
||
| # TODO: r.patch in v.to.rast of values at exact seed site locations. currently set to null | ||
|
|
||
| grass.run_command("r.colors", map=output, color="bcyr", quiet=True) | ||
| grass.run_command( | ||
| "r.support", | ||
| map=output, | ||
| history=" ", | ||
| title="Inverse cost-weighted interpolation", | ||
| ) | ||
| grass.run_command("r.support", map=output, history="v.surf.icw interpolation:") | ||
| grass.run_command( | ||
| "r.support", | ||
| map=output, | ||
| history=" input map=" + pts_input + " attribute column=" + column, | ||
| ) | ||
| # FIXME: if layer !=1 then (layer="$layer") probably needed here | ||
| grass.run_command( | ||
| "r.support", | ||
| map=output, | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
No need to limit the values. Definitively not the upper limit.
1-should work just fine.There was some work done for C and OpenMP (not Python yet) to use zero and negative in a useful way.