[Bug] Discrepancy between outputs of r.what and d.what.rast for the same coordinates

[girishnand]

Describe the bug

I am querying a raster generated from a GeoTIFF file for values at specific lat-long locations. I used d.what.rast and clicked on the displayed raster image using the mouse. The purported value at that location is displayed in a pop-up window, which I redirected to console. I then queried the raster value using r.what, providing it the coordinates returned by d.what.rast (corresponding to the point clicked upon). The raster value returned by r.what is sometimes slightly different from that returned by d.what.rast, even though the coordinates reported are identical.

To Reproduce
Steps to reproduce the behavior:

1. Download and use the GeoTiff file inside the following zip file.
   TL2010_Apr_gf.zip

This is a Linke Turbidity GeoTiff originally downloaded from - https://www.soda-pro.com/help/general-knowledge/linke-turbidity-factor
The precise link for downloading is - https://www.soda-pro.com/documents/10157/326336/Tl2010_MonthlyAv_and_YearlyAv_WGS84_GeoTIFF.zip/b77eb857-29f5-49f7-b7bd-2726f80c54cc

Note: This contains a set of 13 GeoTiff images, out of which I have included the one for April as my test data. I expect the other images to give similar results.

2. Start GRASS GIS and execute the following sequence of commands from the console -

r.in.gdal input=(file-path)\TL2010_Apr_gf.tif output=LinkeApr location=latlong

g.mapset PERMANENT location=latlong

d.mon wx0

d.rast map=LinkeApr

d.what.rast LinkeApr

Click on a random point on the LinkeApr image. When the pop-up window showing co-ordinates and value of the chosen point appears, select the "Redirect to console" checkbox at the lower left corner. When the same output appears in the console, copy the coordinates displayed. run r.what on the console using these self-same coordinates, something like -

r.what LinkeApr coordinates=76.6561513889196,30.37854888375702 -nfir

The following is sample output for two randomly chosen points on the image. The output from d.what.rast (as re-routed to console) appears first, followed by that for r.what for each of the two points. Note that coordinates reported by both commands are identical. -

C:\OSGeo4W\bin>Query results:
east, north: -19.05213269380684, 11.232227483662228
LinkeApr@PERMANENT:
value: 170
color: 170:170:170
--------------------------------------------------

C:\OSGeo4W\bin>r.what LinkeApr coordinates=-19.05213269380684,11.232227483662228
-19.05213269380684|11.232227483662228||170

C:\OSGeo4W\bin>east, north: 76.6561513889196, 30.37854888375702
LinkeApr@PERMANENT:
value: 126
color: 126:126:126
--------------------------------------------------

C:\OSGeo4W\bin>r.what LinkeApr coordinates=76.6561513889196,30.37854888375702 -nfir
easting|northing|site_name|LinkeApr|LinkeApr_label|LinkeApr_color
76.6561513889196|30.37854888375702||136||136:136:136

C:\OSGeo4W\bin>

The value reported for the first point by d.what.rast and r.what is identical, but that for the second point is not.

Expected behavior
d.what.rast and r.what should generate the same output for the same input coordinates.

Screenshots

System description (please complete the following information):

I am on Windows 11 Home & my grass version is 8.3.dev

• details about further software components

C:\Users\Girish>g.version -rge
version=8.3.dev
date=2023
revision=c05cc14387
build_date=2023-04-24
build_platform=x86_64-w64-mingw32
build_off_t_size=8
libgis_revision=a82501dc85
libgis_date=2023-02-27T12:45:25+00:00
proj=9.2.0
gdal=3.6.4
geos=3.11.2
sqlite=3.41.1


C:\Users\Girish>python3 -c "import sys, wx; print(sys.version); print(wx.version())"
3.9.5 (tags/v3.9.5:0a7dcbd, May  3 2021, 17:27:52) [MSC v.1928 64 bit (AMD64)]
4.2.0 msw (phoenix) wxWidgets 3.2.1

[petrasovaa]

The difference is likely caused by resampling, due to using a different computational region, r.what takes computational region into account. When you query in GUI, the region is set based on the display region (meaning, the region the raster was rendered with), that depends on zoom. One reason for GUI to use the display region instead of the current computational region is that you typically want to get the raster value from a place you click at, even though your computational region extent may not include that pixel. But this behavior then leads to the effect you observed, which is not expected either. So perhaps we need to set the region during querying to match the individual rasters (g.region raster=map) and call r.what for each raster.

[girishnand]

@petrasovaa

Using g.region raster=LinkeApr does seem to reduce the frequency of such mis-matches between d.what.rast and r.what, but does not eliminate them entirely. For example -

C:\OSGeo4W\bin>east, north: 67.96208528089304, 27.72511847233084
LinkeApr@PERMANENT:
value: 188
color: 188:188:188
--------------------------------------------------

C:\OSGeo4W\bin>r.what LinkeApr coordinates=67.96208528089304,27.72511847233084 -nfir
easting|northing|site_name|LinkeApr|LinkeApr_label|LinkeApr_color
67.96208528089304|27.72511847233084||184||184:184:184

C:\OSGeo4W\bin>

The first output above is from clicking on the screen after entering d.what.rast LinkeApr on the console, and the second is by using r.what at the console prompt with co-ordinates copy-pasted from the ones reported by d.what.rast. About every fourth or fifth such trial results in mismatched values.

Regardless of region setting, so long as the co-ordinates reported by d.what.rast and r.what are identical, so should the values.

I also loaded the GeoTiff file as an array in GNU Octave and examined the value at the coordinates above (cellsize is 1/12th of a degree along both x and y) -

>>
>> I = imread("TL2010_Apr_gf.tif") ;
>> I( ceil((90-27.72511847233084)*12), ceil( (67.96208528089304+180)*12 ) )
ans = 184
>>

which suggests that r.what is correct, while d.what.rast isn't.

regards,
Girish.

[petrasovaa]

Using g.region raster=LinkeApr does seem to reduce the frequency of such mis-matches between d.what.rast and r.what, but does not eliminate them entirely. For example -

Right, sorry, I was not clear in my comment. It was more of a note to developers, this needs to be fixed in the GUI by setting the region internally.