[GRASS-dev] Re: r.external

[Markus Neteler]

On Sat, Sep 13, 2008 at 2:59 AM, Glynn Clements
<glynn at gclements.plus.com> wrote:
>
> Markus Neteler wrote:
>
>> > r.info -r pat_dtm_5m
>> > WARNING: category support for [pat_dtm_5m] in mapset [PERMANENT] missing
>> > min=nan
>> > max=nan
>>
>> I observe that 'f_range' remains empty:
>>
>> GRASS 6.4.svn (patUTM):~ > l
>> /home/neteler/grassdata/patUTM/PERMANENT/cell_misc/pat_dtm_5m
>> total 52
>> -rw-rw-r-- 1 neteler neteler   132 2008-09-12 14:07 gdal
>> -rw-rw-r-- 1 neteler neteler     0 2008-09-12 14:07 f_range
>> -rw-rw-r-- 1 neteler neteler     5 2008-09-12 14:07 f_quant
>> -rw-rw-r-- 1 neteler neteler    28 2008-09-12 14:07 f_format
>> -rw-rw-r-- 1 neteler neteler 39928 2008-09-12 22:16 histogram
>
> Odd.
>
>> Could the file be populated using GDAL functions?
>
> It is:

(sorry, I was looking only in the library gdal.c)
...

> The fact that the f_quant and f_format files are being written implies
> that G_write_fp_range() is being called, but is writing an empty file.
>
> Hmm; the -r flag won't help if GDAL is reporting a stored range.

I have added debug output:

GRASS 6.4.svn (patUTM):~/grass64/raster/r.external > svn diff
Index: main.c
===================================================================
--- main.c      (revision 33396)
+++ main.c      (working copy)
@@ -303,8 +303,12 @@

     info->range[0] = GDALGetRasterMinimum(hBand, &bGotMin);
     info->range[1] = GDALGetRasterMaximum(hBand, &bGotMax);
+    G_debug(0, "info->range[0]: %f", info->range[0]);
+    G_debug(0, "info->range[1]: %f", info->range[1]);
     if(!(bGotMin && bGotMax))
        GDALComputeRasterMinMax(hBand, !exact_range, info->range);
+    G_debug(0, "exact info->range[0]: %f", info->range[0]);
+    G_debug(0, "exact info->range[1]: %f", info->range[1]);

     G_init_colors(&info->colors);

and it does:

GRASS 6.4.svn (patUTM):~/grass64/raster/r.external > r.external
in=/geodata2/.../pat_DTM_2008_derived_5m_UTM_WGS84.tif out=pat_dtm_5m
--o
Projection of input dataset and current location appear to match
D0/0: info->range[0]: -4294967295.000000
D0/0: info->range[1]: 4294967295.000000
D0/0: exact info->range[0]: 0.000000
D0/0: exact info->range[1]: 3563.149902
<pat_dtm_5m> created

It doesn't match well the GDAL Computed Min/Max (see below) but at
least it is roughly defined:

gdalinfo pat_DTM_2008_derived_5m_UTM_WGS84.tif -mm
...
Computed Min/Max=0.000,3758.871

The file "f_range" continues to be empty.

> But r.support.stats should fix that.
>
> Is it possible that something is causing GDAL to "see" NaNs?
>
> Maybe the presence of NaNs interferes with the range calculation. As
> all inequalities return false for NaN, code such as:
>
>        #define MIN(a,b) ((a) < (b) ? (a) : (b))
>        #define MAX(a,b) ((a) > (b) ? (a) : (b))
>
>        for (i = 0; i < n; i++) {
>            min = MIN(min, x[i]);
>            max = MAX(max, x[i]);
>        }
>
> will set both min and max to NaN if even a single x[i] is NaN.

Adding in create_map():

@@ -454,12 +458,16 @@
        struct Range range;
        range.min = (CELL)info->range[0];
        range.max = (CELL)info->range[1];
+       G_debug(0, "CELL range.min: %d", range.min);
+       G_debug(0, "CELL range.max: %d", range.max);
        G_write_range(output, &range);
     }
     else {
        struct FPRange fprange;
        fprange.min = info->range[0];
        fprange.max = info->range[1];
+       G_debug(0, "FCELL range.min: %f", fprange.min);
+       G_debug(0, "FCELL range.max: %f", fprange.max);
        G_write_fp_range(output, &fprange);
        write_fp_format(output, info);
        write_fp_quant(output);

shows
GRASS 6.4.svn (patUTM):~/grass64/raster/r.external > r.external
in=/geodata2/.../pat_DTM_2008_derived_5m_UTM_WGS84.tif out=pat_dtm_5m
--o
...
D0/0: FCELL range.min: 0.000000
D0/0: FCELL range.max: 3563.149902

So the range (slightly incorrect) survives, but
G_write_fp_range() isn't writing it out.

Markus