(carried here from grass-stats, discussing improvements of the
R-GRASS interface)
On Tue, Jun 24, 2008 at 1:48 PM, Roger Bivand <Roger.Bivand@nhh.no> wrote:
On Tue, 24 Jun 2008, Markus Neteler wrote:
On Tue, Jun 24, 2008 at 12:23 PM, Roger Bivand <Roger.Bivand@nhh.no> wrote:
Is there a shell-level utility in GRASS to check whether the region of a
raster is equal to the current working region - the function could call
that and create a copy of the required raster in current resolution if they
differ?
and
g.region -g | grep '^n=\|^s=\|^e=\|^w='
n=228513
s=214975.5
w=629992.5
e=645012
But the resolution may differ too - I need a one shot comparison. Isn't
anything exposed - I've looked in g.copy and r.resample without seeing
anything. Maybe:
>>> Is there a shell-level utility in GRASS to check whether the region of a
>>> raster is equal to the current working region - the function could call
>>> that and create a copy of the required raster in current resolution if they
>>> differ?
>>
>> The easiest way it to compare the output of
>>
>> r.info -g map127
>> north=230000
>> south=214000
>> east=646000
>> west=628000
>>
>> and
>> g.region -g | grep '^n=\|^s=\|^e=\|^w='
>> n=228513
>> s=214975.5
>> w=629992.5
>> e=645012
>
> But the resolution may differ too - I need a one shot comparison. Isn't
> anything exposed - I've looked in g.copy and r.resample without seeing
> anything. Maybe:
>
> r.info -gs elevation.dem
> g.region -g | grep '^n=\|^s=\|^e=\|^w=\|^ewres=\|^nsres='
>
> but an internal match (where r.resample would be a no-op or not) would be
> helpful.
Any hints for Roger?
"r.info -gs" will print both the bounds and resolution in shell-script
style:
Is there a shell-level utility in GRASS to check whether the region of a
raster is equal to the current working region - the function could call
that and create a copy of the required raster in current resolution if they
differ?
and
g.region -g | grep '^n=\|^s=\|^e=\|^w='
n=228513
s=214975.5
w=629992.5
e=645012
But the resolution may differ too - I need a one shot comparison. Isn't
anything exposed - I've looked in g.copy and r.resample without seeing
anything. Maybe:
Although the output of the two commands is slightly different, it's
close enough that the keys will sort the same, so you can use e.g.:
map=map127
( g.region -g ; r.info -gs $map ) | sort | grep '^[nsew]' | sed 's/^.*=//' | uniq -u | grep . &> /dev/null
if [ $? = 0 ] ; then
echo different
else
echo same
fi
Thanks, elegant, I know. But it isn't cross-platform. All I need is an answer TRUE/FALSE for the comparison of the 6 values. I assume that the printf()s are using the same rounding, right? Both r.info and g.region use G_format_easting(), G_format_northing(), G_format_resolution() in lib/gis/wind_format.c, using format_double(), which is hard-wired to "%.8f". Are we going to pick up numerical fuzz anywhere from conversion to text (high longlat resolution)?
I guess I'll have to trust the text representation, read it into R, convert to doubles, and ask abs(rast$w - wind$w) < .Machine$double.eps for each of the six comparisons, and perhaps handle longlat with a slacker precision.
Roger
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand@nhh.no
>>> but an internal match (where r.resample would be a no-op or not) would be
>>> helpful.
>>
>> Any hints for Roger?
>
> "r.info -gs" will print both the bounds and resolution in shell-script
> style:
>
> $ r.info -gs elevation.dem
> north=4928000
> south=4914020
> east=609000
> west=590010
> nsres=30
> ewres=30
>
> Although the output of the two commands is slightly different, it's
> close enough that the keys will sort the same, so you can use e.g.:
>
> map=map127
> ( g.region -g ; r.info -gs $map ) | sort | grep '^[nsew]' | sed 's/^.*=//' | uniq -u | grep . &> /dev/null
> if [ $? = 0 ] ; then
> echo different
> else
> echo same
> fi
>
Thanks, elegant, I know. But it isn't cross-platform. All I need is an
answer TRUE/FALSE for the comparison of the 6 values. I assume that the
printf()s are using the same rounding, right? Both r.info and g.region use
G_format_easting(), G_format_northing(), G_format_resolution() in
lib/gis/wind_format.c, using format_double(), which is hard-wired to
"%.8f". Are we going to pick up numerical fuzz anywhere from conversion
to text (high longlat resolution)?
Yeah, lat/lon could be a problem (r.info uses d:m:s while g.region
uses decimal), as could fractional coordinates, due to precision
issues.
I guess I'll have to trust the text representation, read it into R,
convert to doubles, and ask abs(rast$w - wind$w) < .Machine$double.eps for
each of the six comparisons, and perhaps handle longlat with a slacker
precision.
If .Machine$double.eps is what I think it is, it's too small. You
would need to use the formatting precision (i.e. 1e-8 for %.8f).
Actually, if you're trying to determine whether r.resample is a no-op,
the ideal test would be to use the region definition to convert the
centres of the corner cells to row/col indices, and check that the
results are (0, 0, cols-1, rows-1). That will catch the case where the
region isn't quite the same as the raster, but close enough that
resampling is actually a no-op.
Yeah, lat/lon could be a problem (r.info uses d:m:s while
g.region uses decimal), as could fractional coordinates, due to
precision issues.
so parse "g.region -p" instead of -g?
I'm making progress with the plugin access to the code underlying gdalinfo for the raster, but getting more precision from g.region for the resolution is problematic. I can however use the rows, columns counts. The longlat raster resolution fuzz is about 0.5e-06 for a region set from a raster and the same value observed through the plugin from the raster and as printed by g.region. Getting closer.
Roger
Hamish
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand@nhh.no