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

[Victor Wren]

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

If the clipping window for the import is determined by the current region 
(the usual Grass behaviour), 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?  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.

> 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.
 
> > m.proj does not seem to support datum conversion (and m.datum.shift doesn't 
> > support UTM coordinates)...
 
> AFAICT, the CVS version of m.proj2 should handle datum conversion.

I fear CVS. :-)  I'll wait for 5.0.2 final.  Interestingly, adding 
"datum=nad??" to the 5.0.1 command string for m.proj2 doesn't result in any 
error, but the results it comes up with disagree with Corpscon by about 46m 
easting (some days I just don't know WHO to believe!).

> > 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).  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'd still love a flag for 
"color outside the lines" so that I can just ignore the current region. :-)  

> > I'm using Grass 5.0.2.
> 
> 5.0.2 hasn't been released yet; it's either 5.0.1, or one of the CVS
> versions (which are currently named "5.0.2-cvs").

Yep, double-checking, it's 5.0.1.  I plead temporal insanity.

Thanks for the help.  Without -n, it still works in an oversized area, but 
the final raster is the correct size, and at least it doesn't try to allocate 
memory for the whole region (which makes the conversion go a LOT faster)!  
It's still not exactly script-friendly, but I can cope, now.

Victor Wren