[GRASS-user] Re: r.in.gdal oddity

[Hamish]

Hamish wrote:
> > WRT the -l flag, if you use 'r.in.gdal -l' you _must_
> > use r.region after to correct the map bounds+resolution.
...
> Maybe I am silly. But use r.region how?

n= s= e= w=

> - with -l, r.in.gdal will import ALL pixels, right?

yes, the actual array data is imported untouched,
just the
north:
south:
east:
west:
in $MAPSET/cellhd/ are changed.

r.in.gdal man page:
"While the resulting bounds and resolution will
 likely be wrong for your map the  map's  data
 will  be  unaltered and safe. After resetting to
 known bounds with r.region you should double check
 them with r.info,  paying  special  attention  to
 the map resolution."


> - then I have them in (even one pixel line/row
> too much perhaps than possible by definition)

which dataset are you looking at?

for the global population count one,
  gl_gpwfe_pcount_10_ascii_one.zip
the cellsize setting in the original dataset was
simply wrong. fix the resolution in the header and
it imports ok.



for the WorldClim precipitation one,
 prec_10m_bil.zip
it is not so clear who's responsible for the problem,
but the thin slice happens because:

$ cat prec2.hdr
BYTEORDER     I
LAYOUT        BIL
NROWS         900
NCOLS         2160
NBANDS        1
NBITS         16
BANDROWBYTES  4320
TOTALROWBYTES 4320
BANDGAPBYTES  0
NODATA        -9999
ULXMAP        -179.9166666666667
ULYMAP        89.91666666666667
XDIM          0.166666666666667
YDIM          0.166666666666667

DatabaseName  WorldClim
Version       1.4
Release         3
Created       20060101
Projection    GEOGRAPHIC
Datum         WGS84
MinX          -180
MaxX          180
MinY          -60
MaxY          90
MinValue      0
MaxValue      736
Variable      Precipitation
month         February
Units         mm


so:
 MinX != ULXMAP   (?!?)

and:
 west: -179.9166666666667 = +180.083333333333
 east: west + NCOLS*XDIM = w+2160*0.166666666666667
             = 180.083333333334

thus due to cumulative rounding error and the way
$east is calculated it sees east>west, and all
columns fitting within that sliver. but that makes
each column coord unresolvable due to the limits of
FP precision and the map unusable until you set
E<=W instead of E>W.

the culprit here seems to be the weird ULXMAP value.
so our question is if that $west value is legitimate or foul?

> - then I correct that with r.region, but in which
> sense? I would need to cut this row or all gets
> messed up?

since in the .bil header we have:
NROWS         900
NCOLS         2160

it is cell-center convention. If it were grid-
confluence convention it would be n+1, 901 x 2161.
(same is true for the population map, row*col is even)

now we also observe the cells are 1/2 cell off from
nice round numbers:
 -179.9166666666667 - 0.166666666666667/2 = -180

so I suspect what happened is that the original
data was based on a grid-confluence convention,
but was exported or cropped to cell-center raster
convention in a funny way before it was distributed.

GRASS65> r.in.gdal -o in=prec1.bil out=WorldClim.prec.january 
WARNING: Over-riding projection check
 100%
r.in.gdal complete. Raster map <WorldClim.prec.january> created.


(no -l flag needed as not exceeding 90deg latitude)

GRASS65> r.info WorldClim.prec.january
...
 |            N:        90N    S:        60S   Res:  0:10                     |
 |            E:       180E    W:       180E   Res:  0:10                     |

at 10' res dataset I see no problem,
 [GRASS 6.5.svn46412 (2011), perhaps slightly old]
but if needed, using:

GRASS65> r.region WorldClim.prec.january e=w+0

should fix it. (e=w doesn't seem to pass the parser)



The data still needs to have its nulls removed,
.bil header:
MinValue      0
MaxValue      885

r.info:
minimum: 0
maximum: 55537

GRASS65> r.null WorldClim.prec.january setnull=55537

GRASS65> g.region rast=WorldClim.prec.january
GRASS65> r.colors -e WorldClim.prec.january color=bcyr
GRASS65> r.univar WorldClim.prec.january
...
maximum: 736


but looks ok to me..


regards,
Hamish