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

[Glynn Clements]

Victor Wren wrote:

> > > Given that I won't know the extents of a raster in the current map until I 
> > > import it, it seems a chicken-and-egg problem to decide what the current 
> > > region should be set to before importing. 
> 
> > This shouldn't be the case; r.proj includes code to project the
> > borders of both the current region and the source raster to determine
> > the final borders.
> 
> According to the man pages, the output limits of the map are determined by 
> the current region.  Does that mean that the current region is supposed to be 
> a clipping function (i.e. the map will be bounded by the current region, but 
> not expanded to fit it?)  This is still not very friendly where scripting is 
> concerned (especially given that SDTS for DEM is generally NAD27 and SDTS for 
> DLG data is more often NAD83 [maybe] meaning something almost always has to 
> be project to combine them)  So say I want run a script to project DEM 
> rasters scattered all over the  state of Colorado from nad27 to nad83 at 
> 1:24,000 -- should I set my region to cover the whole state (plus about 
> 20,000m for outliers) before I start?

Probably. Assuming that r.proj works correctly, the resulting region
shouldn't be any larger than necessary. OTOH, if the current region is
initially too small, it won't be enlarged (to allow you to extract a
section of a larger map).

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

> If the clipping window for the import is determined by the current region 
> (the usual Grass behaviour),

Import utilities (r.in.*) ignore the current region. r.proj may shrink
the current region to reduce the processing time.

> I would need to have a close guess for each of 
> my area DEMs before I can bring them over to keep them from getting trimmed, 
> or pick a region that I know will exceed the dimensions of any DEM I am 
> projecting.  Yes?

Yes.

> If region is that vital to import, I would suggest there 
> should be SOME way to glean the bounds needed to set the region from mapsets 
> in other locations, without having to hand-carry it (something like "g.region 
> rast= zoom=" combined with r.proj)  It would be cool if this could be 
> automagic, just by adding a mapset, database and location parameters to 
> g.region.

It shouldn't be that hard to add an option to r.proj to enlarge the
destination region to encompass the entire source map.

> > What should happen is that r.proj projects the current region to the
> > source location's projection to determine the region of the source map
> > which needs to be read, then projects that region back to the
> > destination location to determine the extents of the resulting map.
> 
> That's an interesting approach. I presume this is to theoretically
> minimize the amount of data that has to be converted.

It's done for both speed and memory.

The relevant part of the source map has to be kept in memory, so it's
preferable not to read portions which won't actually be used. The
projection is performed by iterating over all cells in the (final)
destination region, so it's preferable to shrink the destination
region as small as possible.

There are a few (extreme) cases where region calculation falls down
(e.g. polar regions), so the -n switch exists to disable it (r.proj
will read the entire source map into memory, and iterate over the
entire destination region).

> > > A switch in [rsv].proj to override the current region and create a new raster 
> > > with the same geometric extents (projected) as the source raster would be a 
> > > real timesaver, but I'm not enough of a code hog to even know whether that 
> > > can be done at the level that r.proj works at.
> > 
> > Provided that the current region is initially larger than the
> > resulting map (and you don't use the -n flag), you should end up with
> > a map whose extents are "shrink wrapped" to the actual data. If that
> > doesn't happen, it's a bug.
> 
> Ahah, Eureka, gadzooks, et al!  The "-n" flag was the problem.  The flag 
> doesn't exist in the man page, but in r.proj -help, it says "don't attempt to 
> crop output map" which sounded like what I wanted.  That seems a little 
> misleading to me, because even with the -n flag, the output is STILL cropped 
> to the current region.  All that -n seems to do is use the current region for 
> the output raster dimensions (whether it's bigger OR smaller than the source 
> raster).

Yes.

BTW, "-n" shouldn't make the results any worse (i.e. you shouldn't
lose even more data; it may make things worse in terms of CPU, memory
and disk usage, though).

> 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. That's rather verbose, and
unclear (particularly if you haven't got used to how other GRASS
programs work); Unless someone comes up with something better, I'll
probably settle for "disable optimization".

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.

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

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