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

[Victor Wren]

> > 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 shouldn't be that hard to add an option to r.proj to enlarge the
> destination region to encompass the entire source map.

I would hope that's true.  There's an option like that in r.in.gdal already.

> > 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.

> programs work); Unless someone comes up with something better, I'll
> probably settle for "disable optimization".

I could go for such a terse explaination, but only if there is an explanation 
in the man page what that amounts to.  I'd be happy to write it, but I don't 
know the generalized ramifications, since I've been working in a very limited 
set of projections.

> Basically, the "shrink wrapping" code works by projecting the
> *borders* of the region in question, returning the bounding box of the
> projected points. This falls down if you e.g. project a polar region
> to a cylindrical projection, so I added an option to turn it off.
> 
> > I'd still love a flag for 
> > "color outside the lines" so that I can just ignore the current region. :-)  
> 
> Existing behaviour is consistent with the model that all maps have
> infinite extent (even if everything outside of an arbitrary finite
> rectangle happens to be "no data"), and it's up to the user as to
> which finite section they're interested in.

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.in.??? 
doesn't do region trimming (at least from what I can see in the man pages).  
operational functions like r.mapcalc  are trimmed, but they don't alter 
coordinates (only cell values).  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? 
Or is r.slope.aspect an oddball that defies the GRASS philosophy that nobody 
but me really uses?

> 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.  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)

Now, back to figuring out what an SDTS "manifold" is.  I wonder if it could 
help my horsepower...

Victor Wren





Victor Wren
Designer,
Timension Inc.
1350 C Pear Ave
Mountain View CA 94043
(650) 564-9397
Fax: (650) 564-9398
Opinions stated in this letter are not necessarily
those of Timension Inc. or the management.  All
Rights Reserved.  No spitting.