[GRASSLIST:5854] Question regarding r.proj

I need a suggestion how to proceed. I've got a large number of raster maps
that I need to project from NAD27 datum to NAD83 (both in UTM, grs80
ellipsoid, though I think NAD27 is supposed to be clark66, but corpscon
doesn't support that ellipsoid, so I couldn't check the coordinate shift
there). I wrote a shell script to do this, but it failed dismally, because I
didn't have the correct region selected for each raster at the time of
import.

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. If the region is too large, the
raster will be expanded on import to its limits (requiring a huge amount of
memory and conversion time). If it's too small, the maps will end up being
trimmed (or will be skipped entirely, if they're outside the region.) I
could set them huge, and r.patch them after, but that seems inelegant (and
I'm not sure how to automate that, either).

m.proj does not seem to support datum conversion (and m.datum.shift doesn't
support UTM coordinates), so I can't readily use that to find out where the
projected raster extents will end up. Experiments with corpscon show that in
my area of interest, the raster will shift about 200m Northing and about 50m
Easting.

I'd really like to automate the movement of rasters between my nad27 and
nad83 maps and do it all inside Grass, if possible. I have a lot of them, and
going through manually, finding the NSEW extents for each of them in nad27,
taking them to corpscon and projecting them to nad83, then setting each
region manually in Grass before doing r.proj is going to be many hours more
work than I was hoping for.

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.

Suggestions very welcome.

I'm using Grass 5.0.2.

Victor Wren

Victor Wren wrote:

I need a suggestion how to proceed. I've got a large number of raster maps
that I need to project from NAD27 datum to NAD83 (both in UTM, grs80
ellipsoid, though I think NAD27 is supposed to be clark66, but corpscon
doesn't support that ellipsoid, so I couldn't check the coordinate shift
there). I wrote a shell script to do this, but it failed dismally, because I
didn't have the correct region selected for each raster at the time of
import.

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. If the region is too large, the
raster will be expanded on import to its limits (requiring a huge amount of
memory and conversion time). If it's too small, the maps will end up being
trimmed (or will be skipped entirely, if they're outside the region.) I
could set them huge, and r.patch them after, but that seems inelegant (and
I'm not sure how to automate that, either).

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.

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.

m.proj does not seem to support datum conversion (and m.datum.shift doesn't
support UTM coordinates), so I can't readily use that to find out where the
projected raster extents will end up. Experiments with corpscon show that in
my area of interest, the raster will shift about 200m Northing and about 50m
Easting.

AFAICT, the CVS version of m.proj2 should handle datum conversion.

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.

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

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

On Tue, 25 Mar 2003, Glynn Clements wrote:

AFAICT, the CVS version of m.proj2 should handle datum conversion.

Yes, but only the cmd version. 'datum=xxx' should be specified as one of
the command-line arguments for both inproj and outproj for datum
transformation to take place. The interactive version does not prompt for
datum, but the code is such a mess (clone of g.setproj) that it is hardly
worth changing without re-writing it totally. Better to use an external
program like cs2cs from the PROJ.4 distribution to transform point
co-ordinates.

Paul

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

> 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

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

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@virgin.net>

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

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

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. :slight_smile: 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@virgin.net>

I've been playing around with r.proj for a couple of days, now. I've noticed
some things that explain unusual behavior that I was seeing in projected DEM
maps. I was noticing that they didn't meet up at the edges quite right after
projection. I figured out what was going on when I decided to see if I could
find a simple way to shrink maps down if the -n flag is necessary. Using
g.region rast={map} zoom={map} was not yielding a map with dimensions equal
to the source map (being that both maps were UTM, they should have been). It
seemed to grow a cell or two in all directions.

What I think is happening is that using any interpolation method will cause
boundary effects, where the edges are being interpolated with nulls
(interpreted as zero?). This has the effect of creating cells where there
were none before (if the region is larger than the source map, and the -n
switch is used), but more importantly, it causes distortions at the edge of
any maps, even when the -n switch isn't used. I've decided for my
application, bilinear interpolation was unnecessary, but in a general case it
might not be (such as projecting from geo coordinates to UTM)

Is there a way to overcome this? I'm a mathematical idiot, but it seems to
me that if any of the cell values in the interpolation group is null (the
four cells of the bilinear group, or the 16 cells of the cubic group), the
interpolation should not be done, or perhaps should be "faked" by repeating
the edge cells out until a legitimate group can be created, which may be no
more valid, but at least will keep the edge from dropping away). If the
center cell is a null, does it get filled in by interpolated data? My
experience would suggest that it does (that would also distort data in
adjacent that is depending on that null cell as one of its bilinear or cubic
group). Should "no-data" be treated as zero for the purpose of
interpolation? I can think of arguments both ways. Does this matter to
anybody?

Victor Wren

I'm trying to bring in a large number of elevation files, in preparation for
turning them into a surface map. The file is about 150,000 floating-point
locations in this format:
2887690.49718 1673263.622019 10969.170819 @" "
2887694.185556 1673257.944703 10969.240661 @" "
2887640.301113 1673223.419907 10972.161139 @" "
.......

I'm importing them to a sites file with:
s.in.ascii input={infilename} sites={sitesname} d=3

This goes well, and when I do d.sites, they show up where they're supposed to
be.

I then want to convert them to a raster map. This is the command I use:

s.to.rast input={sitesname} output={rastname} field=dim findex=3

(I've tried without field and findex, and it doesn't like it I'm not quite
sure exactly what I'm doing, but it seems to work.) It says:

Using size option 0
Sites map Type: 0, Dims: 3, Strs: 1, Dbls: 0

It warns me about dropping multiple sites on one cell (which is fine -- the
sites are overly dense, and the elevations should be very close), and
crunches along for awhile, then segfaults.

I thought it might have a problem with the number of sites, so I trimmed it
back to about 48,000 and it still segfaulted. I trimmed it again 1,000, and
it successfully finished. I only have to do this 600 more times, and I'll be
done with the first step.

What's the practical limit to the number of sites that s.to.rast can cope
with?

Victor Wren

On Wed, 26 Mar 2003, Glynn Clements wrote:

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.

should note here that setting the resolution with the command line option
resolution= is more or less obsolete. a better way of getting _exactly_
the resolution you want is to create a region within your location with
the desired resolution, make it the current region and then project into
that. (of course, you can't get anything better than, not even equal to,
the input map. a raster projection, contrary to a vector projection always
involves resampling)

another thing worth pointing out in this thread is that the only usable
interpolation method for DEM maps is 'nearest'. both bilinear and cubic
give severe distortions of coastlines, shape of lakes and other outlines
around flat surfaces.

Morten Hulden

Victor Wren wrote:

I've been playing around with r.proj for a couple of days, now. I've noticed
some things that explain unusual behavior that I was seeing in projected DEM
maps. I was noticing that they didn't meet up at the edges quite right after
projection. I figured out what was going on when I decided to see if I could
find a simple way to shrink maps down if the -n flag is necessary. Using
g.region rast={map} zoom={map} was not yielding a map with dimensions equal
to the source map (being that both maps were UTM, they should have been). It
seemed to grow a cell or two in all directions.

r.proj explicitly grows the region by two cells in each direction, to
allow for interpolation.

What I think is happening is that using any interpolation method will cause
boundary effects, where the edges are being interpolated with nulls
(interpreted as zero?). This has the effect of creating cells where there
were none before (if the region is larger than the source map, and the -n
switch is used), but more importantly, it causes distortions at the edge of
any maps, even when the -n switch isn't used. I've decided for my
application, bilinear interpolation was unnecessary, but in a general case it
might not be (such as projecting from geo coordinates to UTM)

Is there a way to overcome this? I'm a mathematical idiot, but it seems to
me that if any of the cell values in the interpolation group is null (the
four cells of the bilinear group, or the 16 cells of the cubic group), the
interpolation should not be done, or perhaps should be "faked" by repeating
the edge cells out until a legitimate group can be created, which may be no
more valid, but at least will keep the edge from dropping away). If the
center cell is a null, does it get filled in by interpolated data? My
experience would suggest that it does (that would also distort data in
adjacent that is depending on that null cell as one of its bilinear or cubic
group). Should "no-data" be treated as zero for the purpose of
interpolation? I can think of arguments both ways. Does this matter to
anybody?

AFAICT, the cubic interpolation code returns null if the centre cell
is null, and replaces null cells by the contents of the centre cell.

This looks like a fudge; personally, I think that the result should be
null if any of the cells used by the interpolation are null. (Either
that, or there should be separate interpolation algorithms for each
possible combination of null/non-null cells).

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