[GRASSLIST:5882] Re: Question regarding r.proj

[Glynn Clements]

Victor Wren wrote:

> > > Since output resolution can be 
> > > specified by r.proj, could I set the region resolution to be something like 
> > > 5000m before importing?  (yah, yah, I know -- try it and see)
> > 
> > Only the final resolution (set by the current region in the
> > destination location or the resolution= option to r.proj) should
> > matter. Import utilities import cell-for-cell (i.e. without any
> > rescaling), and r.proj reads the source map cell-for-cell.
> 
> In practice, this is almost true.  I've experimented with it, and the 
> resolution is only correct IFF the resolution of the destination region is 
> the same as the source map.  For chuckles, I set the resolution of my region 
> to ns-15m and ew-56m (it did not come out even when it recalculated to make 
> the intervals fit the bounds, of course, since I pulled the numbers out of a 
> hat).  When I project a DEM10 nad27 raster into this less-than-perfect--nad83-
> world, r.proj reports the source resolution as 10 ew and 10 ns, but the 
> destination resolution is 9.999863 ew and 10.003091 ns.  If I reset the 
> region to be a nice, even 10x10 resolution, the destination resolution is 
> also 10x10.  Bug?

It looks like it.

> > > In other words, a more accurate description for the -n flag would 
> > > be "Use current region for raster dimensions in preference to source raster 
> > > dimensions, padding with nulls if necessary."
> > 
> > I suppose "don't crop output map any more than any other GRASS program
> > would do" is a more accurate description.
> 
> Is that accurate?  I have not encountered another GRASS program that will 
> expand the output to create a raster that encompasses the entire current 
> region.  I'm aware that a smaller (or off-centered) region will crop data 
> with most operations (this usually bites me with r.mapcalc), but I wasn't 
> aware that dealing with a large region by adding null padding cells is SOP.  
> I'm still new at this, of course, and I should play around to see if I know 
> what I'm talking about.

Yes, it's accurate.

When a program creates a new map, it will have the bounds and
resolution taken from the region. When a program opens an existing
map, it is automatically cropped, padded and rescaled according to the
region settings.

This is true for any program which doesn't explicitly change the
region, so it's true for most raster programs.

> I'm interested in the section I'm importing. :-)  Seriously, though, r.proj 
> functionally appears to me to be in the family of import functions.

r.proj definitely isn't an import module; it's more like a
conventional processing module, except that the source map is
typically in a different location.

> r.in.??? doesn't do region trimming (at least from what I can see in
> the man pages).

The r.in.* programs don't generally follow the normal behaviour
(crop/pad/rescale to the current region). Instead, they normally
import cell-for-cell, and set the imported map's region so that one
corner is at the origin of the coordinate system (this way, it should
be obvious that its bounds are wrong). You then have to use r.region
or r.support to set the correct parameters.

The exception is where georeferencing information is present, e.g. a
.tfw file, or GeoTIFF; in that case, the import utility should set the
map's bounds and resolution correctly.

> operational functions like r.mapcalc  are trimmed, but they don't alter 
> coordinates (only cell values).

The bounds and resolution of the output map will match the current
region. Any input maps will be cropped, padded and/or rescaled
according to the current region. All processing reads from and writes
to identically sized and aligned grids which match the current region.

> I note in r.slope.aspect that it says 
> 
> "To ensure that the raster elevation map layer is not inappropriately 
> resampled, the settings for the current region are modified...
> If the user really wants the elevation map resampled to the current region 
> resolution, the -a flag should be specified."
> 
> Uh.... How consistent is this region cropping behavior across GRASS
> programs? 

It isn't.

> Or is r.slope.aspect an oddball that defies the GRASS philosophy

Yes; although it certainly isn't the only one.

> that nobody but me really uses?

I suspect that r.slope.apect is commonly used, but it's default
behaviour (without "-a") isn't "normal" GRASS behaviour. It's
behaviour with "-a" is normal GRASS behaviour.

> > The "boundaries" of a map (according to e.g. r.info) are mostly an
> > implementation issue (i.e. disk space). For consistency, we should
> > probably get rid of "g.region rast=..." (OK; maybe not).
> 
> You do and I'll scream.  I use it enough I need a hotkey for it.  "g.region 
> rast=" is a beautiful example of how things can be made to make life easier, 
> by taking data from one place, and putting it in another place, without 
> making GIS weenies like me have to keep a notebook and calculator handy.  I 
> may be missing something, here.

Note that you can have named regions, i.e. "g.region region=..." and
"g.region save=...".

In the general case, the boundaries of your data sources and your
region of interest are two different things. The data sources may
cover a far larger region than you are interested in, or you may have
to patch multiple sources together to cover your region. The only
thing that's really certain is that any given source will overlap your
chosen region to some extent (otherwise you wouldn't bother using it).

> Are you saying that there's not actually a 
> grid somewhere filled in with NULLs when the boundaries are expanded? (and 
> how deep did I want to get into file structures... Hmm)

A GRASS raster is stored as a simple grid. However, the functions
which are used to read a map automatically resample the map according
to the current region settings. So, conceptually it behaves more like
a bivariate function f(x,y) over R^2 than a discrete grid.

If a program has a reason to read the raster cell-for-cell, without
cropping, padding or rescaling, it has to explicitly set the region to
match the raster. Some programs may do this (e.g. r.slope.aspect), but
it isn't the "normal" behaviour.

-- 
Glynn Clements <glynn.clements at virgin.net>