[Bug] Histogram in d.legend does not correspond to actual range of values inside smaller region

[NikosAlexandris]

Describe the bug The histogram added via -d using d.legend, on top of the legend of a raster map, does not correspond to the raster values of the map in the current region. d.histogram, however, correctly plots the distribution of values of the same raster map in the current region. I am unsure if this is a proble with d.legend or related to the NULL file of the raster map.

Expected behavior

d.legend -d should create a histogram that matched the range of the distribution of a raster map's values inside the current computatinal region.

To Reproduce In the following Location

→ g.proj -jf
+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +no_defs +a=6371007.181 +b=6371007.181 +nadgrids=@null +type=crs  +to_meter=1

set the region

→ g.region raster=MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled -gf
projection=99
zone=0
n=5559752.598333 s=4447802.078667 w=0 e=1111950.519667 nsres=926.62543306 ewres=926.62543306 rows=1200 cols=1200 cells=1440000

and plot the legend of a MODIS product

→ d.erase && export GRASS_RENDER_HEIGHT=350 && d.legend MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled title=MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled at=25,35,5,95 fontsize=25 -s -d && tygrass
Drawing horizontal legend as box width exceeds height

So far good. Next, zoom in 'Freiburg'

→ g.region vector=freiburg -gf
projection=99
zone=0
n=5346570.04676644 s=5326179.80211189 w=567744.74363547 e=594970.2834827 nsres=926.82930248 ewres=938.81171887 rows=22 cols=29 cells=638

and re-plot

→ d.erase && export GRASS_RENDER_HEIGHT=350 && d.legend MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled title=MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled at=25,35,5,95 fontsize=25 -s -d && tygrass
Drawing horizontal legend as box width exceeds height

The histogram is obviously (when looking at the map itself) wrong. Some meta:

→ r.describe -d MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled
 100%
* 0.955514-0.960643 0.962353-0.972612 0.974322-0.984580 0.986290-0.988000

→ r.univar MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled -g
n=633
null_cells=5
cells=638
min=0.958
max=0.988
range=0.03
mean=0.979977883096371
mean_of_abs=0.979977883096371
stddev=0.00458923982902483
variance=2.10611222083078e-05
coeff_var=0.468300347200134
sum=620.326000000003

and the map itself with values and the legend from the previous d.legend command:

Tracking back the map's history:

→ r.info -h MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled
Data Source:

Data Description:
   generated by r.mapcalc
Comments:
   MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31 * 0.002 + 0.49

and the original map

→ r.info -h MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31
Data Source:

Data Description:
   generated by r.in.gdal
Comments:
   r.in.gdal input="HDF4_EOS:EOS_GRID:"MYD21A1D.A2019158.h18v04.006.201\
   9202063141.hdf":MODIS_Grid_Daily_1km_LST21:Emis_31" output="MYD21A1D\
   .A2019158.h18v04.006.2019202063141.Emis_31" memory=300 offset=0 num_\
   digits=0

Obviously, it is not an external (r.external) nor a reclassified (r.reclass) map. Before the scaling applied via r.mapcalc on the original map, I did re-set 0 values to become NULL (in a for loop for all maps of the series, including the raster map in question)

# glr is an alias for `g.list raster`
for MAP in $(glr pattern=MYD21A1*Emis_??) ;do echo r.null map=$MAP setnull=0 ;done

Create/reset the null file for this raster map

→ r.support MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled -n
Writing new null file for
[MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled]...
 100%

→ r.info -r MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled
min=0.552
max=0.988

and updating the histogram/range

→ r.support MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled -s

Updating statistics for
[MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled]
Updating the number of categories for
[MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled]

→ r.info -r MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled
min=0
max=1

This seems to have fixed the histogram, yet the range isn't the one we'd expect

0 is added back:

→ r.describe -d MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled
 100%
0.000000-0.003922 0.952941-0.984314

→ r.univar MYD21A1D.A2019158.h18v04.006.2019202063141.Emis_31_Scaled -g
n=638
null_cells=0
cells=638
min=0
max=0.988
range=0.988
mean=0.972297805642638
mean_of_abs=0.972297805642638
stddev=0.0865344932487735
variance=0.00748821852182202
coeff_var=8.89999882202539
sum=620.326000000003

Anyhow, playing with r.support and r.null I seem to come back to the same legend-problem when zoomed in freiburg.

I wonder if d.legend -d does not read the range of raster values of a map in the current region. Is the function Rast_read_fp_range() the one in question here, used here https://github.com/OSGeo/grass/blob/53eda832018485b0d02f94755c8cca9c499c528d/display/d.legend/histogram.c#L51?

r.describe with its -d flag, does https://github.com/OSGeo/grass/blob/53eda832018485b0d02f94755c8cca9c499c528d/raster/r.describe/describe.c#L42-L44

System description (please complete the following information):

• GRASS GIS version

  → g.version -rb
  GRASS 7.8.5dev (2020)
  
  ./configure  --with-blas --with-lapack --with-cxx --with-includes=/usr/include/ --with-libs=/usr/lib64/ --with-proj --with-proj-includes=/usr/local/include/ --with-proj-libs=/usr/local/lib/ --with-proj-share=/usr/local/share/proj/ --with-geos --with-geos=/usr/bin/geos-config --with-gdal=/usr/local/bin/gdal-config --with-x --with-motif --with-cairo --with-opengl-libs=/usr/include/GL --without-ffmpeg --with-python=yes --with-wxwidgets --with-freetype=yes --with-freetype-includes=/usr/include/freetype2/ --with-odbc=yes --with-sqlite=yes --with-postgres=yes --with-postgresql=yes --with-postgres-includes=/usr/include/postgresql --with-opencl --with-opencl-includes=/usr/include/CL/ --with-openmp --with-pthread --with-lapack --with-fftw --with-readline --with-regex --with-nls --with-jpeg --with-tiff --with-png --with-netcdf --without-opendwg --enable-largefile=yes --with-bzlib --with-bzlib-includes=/usr/include/ --with-zstd
  libgis revision: 744fcaefa
  libgis date: 2020-10-16T16:00:53+00:00

[HamishB]

Hi Nikos,

d.legend's -d flag is based on the code from d.histogram, so in theory both should give the same output. Moreover both run 'r.stats' in the background to generate the step class data. Does r.stats give the expected output?

My intention was that the -d flag should respect the current computational region.

[NikosAlexandris]

( Hamish, nice to read you 😄 . I've been windmilled with other works (no GRASS GIS at the moment :-/), so I don't even remember the details anymore. Hopefully someone else can pick this up and help. Cheers! )

[HamishB]

Nice to hear from you too. :-) I can have a look at this at some point, you left enough details to follow it up. One thing I'm currently trying to debug is some weirdness around the region resetting itself as part of the display module internal run-twice design. Which is to say if r.stats output looks ok but d.histogram &/or d.legend -d don't, a good check will be to see if this works properly in GRASS 6 but not in GRASS 7+. If so then we can look for a fix in the display mechanism. Might also test direct rendering vs. d.mon output to see if that makes a difference.