r.proj - new script for re-projecting raster-maps

[kang at risca.rs.msu.edu]

Simon,
Sorry, I probably didn't make it clear about my procedure.  My procedure is 
actually from ERDAS.  I remember in one of the ERDAS command that allows one
to project raster map between different coordinate system.  It first ask user
to lay out a grid.  The intersection points of the grid (a site file) will be
used as the POINTS to build the transformation equation.  However, I have to
convert the site file to vector file since we only have "v.proj".  So, once
these POINTS were projected to the target coordinate system, I have two sets
of POINTS in different coordinate systems.  Their coordinates is the body of th
e
POINTS file.

I agreed that i.rectify may be wrong but not really very far off.  It would be
interesting to use my method to project back and compare the differences.
As you say, you still use a re-sample procedure in your script, which may not b
e
better than the rectification - since the grid is even distributed across the
entire region and you can do 2nd or 3rd order with i.rectify2.  If we just talk
about points, i.rectify definitely is wrong.  However, we are talking about 
a two dimensional object - a cell.  Let's say there is a one row map with
512 columns.  The first cell of the row is (x1,y1) and the last cell is (x2,y2)
.
So, the distance between the first and last cell is x=x2-x1, and y1 should be
the same as y2.  Now, you project the two points to another coordinate system.
The new coordinates for the first and last cell are (nx1,ny1) and (nx2,ny1).
OK, will the distance between nx1 and nx2 be the same?  We all know that the
answer is no, nx2-nx1 will not equal to x2-x1.  So, you have to resample the
cells - either change the cell resolution or decrease/increase the number of 
cells in the row.

It would be interesting to hear from the experts to comment on this - what is
the "right" procedure for raster map projection?

Kang,
Michigan State University