[GRASS-dev] [GRASS GIS] #1424: Profile tool always casts to integers (incl. csv output)

GDAL · August 20, 2011, 6:38am

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------
Hi,

as reported on the grass-user mailing list yesterday, export to .csv
always exports as integers, regardless of if the source data is CELL,
FCELL, or DCELL.

Moritz prepared a patch:
{{{
--- profile.py 2011-08-19 13:28:17.000000000 +0200
+++ profile_new.py 2011-08-19 13:28:43.000000000 +0200
@@ -598,7 +598,7 @@
                                    caption=_("Error"), style=wx.OK |
wx.ICON_ERROR | wx.CENTRE)
                      return False
                  for datapair in r['datalist']:
- file.write('%d,%d\n' %
(float(datapair[0]),float(datapair[1])))
+ file.write('%f,%f\n' %
(float(datapair[0]),float(datapair[1])))

file.close()
}}}

but I think we need to go deeper than that, as `%d` shows up all around
the script, both in preparing the source points to feed into r.profile,
and in processing the results of it.

* Rather importantly in lat/lon locations casting the tie points to int
can silently return the wrong results, with coords shifted a long way from
where you expected.

On the input side it will be easy enough to change `%d` to `%.15g`
(technically `%.9g` would be enough to cover sub-millimeter in lat/lon)
for r.profile and r.what's coord input.

On the output site the distance along transect should be at least to the
mm (`%.3f` or `%.3g` assuming the tool preserves r.profile's column 1
distance as meters), and test for map type, exporting column 2 as `%d` for
CELL, `%.7g` for FCELL, and `%.15g` for DCELL. (see #335)

another thing to keep in mind is that map_units is not always going to be
feet, meters, or degrees. there's little reason it couldn't be microns or
parsecs, so to avoid others being limited by our humble imaginations we
should avoid unnecessary assumptions about that whenever possible.

Hamish

ps- see also #1292 "Profile tool incorrectly processes no data values"

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424>
GRASS GIS <http://grass.osgeo.org>

GDAL · August 20, 2011, 7:03am

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------

Comment(by hamish):

Also I notice an occurrence of `sqrt(a^2 + b^2)` in profile.py, which is
to be avoided as it gives the wrong answer from lat/long locations.

{{{
if lasteast and lastnorth:
dist = math.sqrt(math.pow((lasteast-point[0]),2) + math.pow
((lastnorth-point[1]),2))
}}}

in trunk there is the m.measure C module to use a front end to
G_distance().

In 6.x we used to have swig/python/examples/m.distance as a python front-
end to that; I'm not sure if a simple ctypes replacement exists yet or
not.

thanks,
Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424#comment:1>
GRASS GIS <http://grass.osgeo.org>

GDAL · August 20, 2011, 7:05am

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------

Comment(by hamish):

math.sqrt() is also used in `mapdisp_window.py` and `nviz_mapdisp.py`, and
will need to be checked/replaced.

Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424#comment:2>
GRASS GIS <http://grass.osgeo.org>

GDAL · August 20, 2011, 9:42am

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------

Comment(by glynn):

Replying to [comment:1 hamish]:
> Also I notice an occurrence of `sqrt(a^2 + b^2)` in profile.py, which is
to be avoided as it gives the wrong answer from lat/long locations.

That will give the wrong answer for all locations; Python follows C and
uses `^` for the bitwise exclusive-or operator. Exponentiation uses `**`
(which should be preferred over math.pow()), i.e. sqrt(a**2 + b**2).

However, Python's math module also has a hypot() function, which is
probably the optimal approach here.

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424#comment:3>
GRASS GIS <http://grass.osgeo.org>

GDAL · August 20, 2011, 1:16pm

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------

Comment(by hamish):

Replying to [comment:3 glynn]:
> That will give the wrong answer for all locations; Python follows C and
> uses `^` for the bitwise exclusive-or operator. Exponentiation uses `**`
> (which should be preferred over math.pow()), i.e. sqrt(a**2 + b**2).

oh, I was just pseudo-coding. The actual python code uses math.pow(,2)

> However, Python's math module also has a hypot() function, which
> is probably the optimal approach here.

but not if the map units are in degrees, which is why I'm looking for a
python front-end to G_distance(), or failing that a pipe_command() to the
m.measure module.

finishing a patch for the `%d` problem,
Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424#comment:4>
GRASS GIS <http://grass.osgeo.org>

GDAL · September 23, 2011, 4:01am

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------

Comment(by hamish):

basic patch added to the ticket. But before applying in svn:
`grass.raster_info(raster)['datatype']`
should be used to choose export precision based on the raster map's type,
something like:

{{{
@@ -585,9 +585,18 @@
                                    message = _("Unable to open file <%s>
for writing.") % pfile,
                                    caption = _("Error"), style = wx.OK |
wx.ICON_ERROR | wx.CENTRE)
                      return False
+
+ rast_type = grass.raster_info(raster)['datatype']
+ if rast_type == 'CELL':
+ write_fmt = '%.3f,%d\n'
+ if rast_type == 'FCELL':
+ write_fmt = '%.3f,%.7g\n'
+ if rast_type == 'DCELL':
+ write_fmt = '%.3f,%.15g\n'
+
                  for datapair in r['datalist']:
- file.write('%d,%d\n' %
(float(datapair[0]),float(datapair[1])))
-
+ file.write(write_fmt %
(float(datapair[0]),float(datapair[1])))
+
                  file.close()

dlg.Destroy()
}}}
(vs trunk; untested)

... but why are we doing that at all? for the .csv export why not just re-
run r.profile and save to the output= file directly, with all the above
stuff already taken care of? no risk of introducing aliasing or bugs if we
produce that straight from the source + it would make the code a bit
simpler.

a good test for this is to go into Spearfish, and prepare a limited range
map:
{{{
g.region rast=elevation.dem
r.mapcalc "elev_1k = elevation.dem / 1000.0"
r.info -r elev_1k
min=1.066
max=1.84
}}}

then pull that into a lat/long location with r.proj. So both map units
(decimal degrees) and z-values will show up as obviously wrong if
something gets cast to an integer.

Hamish

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424#comment:5>
GRASS GIS <http://grass.osgeo.org>

GDAL · September 6, 2013, 10:10am

#1424: Profile tool always casts to integers (incl. csv output)
-------------------------------------+--------------------------------------
Reporter: hamish | Owner: grass-dev@…
Type: defect | Status: new
Priority: major | Milestone: 6.4.2
Component: wxGUI | Version: svn-releasebranch64
Keywords: profile tool, precision | Platform: All
Cpu: All |
-------------------------------------+--------------------------------------

Comment(by mlennert):

A student of mine was bitten by this again, recently, so I looked into it.

I'll attach a new patch which combines Hamish' original patch with the
addition of using m.measure. This solves the issue of casts to integers
(and no need for all the raster format checks Hamish suggested as the data
actually '''is''' the original r.profile output with just nulls filtered
out. And .15g seems appropriate for me as the output is then adapted to
the type of data in input (i.e. integers are displayed as integers).

I'm still a bit stuck on the use of the profile tool in lat-long
locations. The use of m.measure was only one aspect. Antoher issue is on
line 232ff where the resolution to be fed to r.profile is calculated.
Using m.measure to determine the transect length gives a length in meters
and thus the calculated resolution is in meters which then does not work
for r.profile in a lat-long location.

So, the patch is still unfinished work, but I'm very busy with other
things these days, and do not have the time in the immediate future to
look at this. Maybe someone else can pick up on what I did.

In any case, it would be great to have the integer->float parts of the
patch applied (and backported to 64release ?).

Moritz

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/1424#comment:6>
GRASS GIS <http://grass.osgeo.org>