[Gdal-dev] Re: [GRASSLIST:6144] r.proj and gdawarp give differernt results

[Frank Warmerdam]

On Sun, 13 Mar 2005 19:41:30 +0000 (GMT), Paul Kelly
<paul-grass at stjohnspoint.co.uk> wrote:
> r.proj is more a backward projection where you have a region of interest
> in a new co-ordinate system and you want to take the portion of the image
> that overlaps with your new system and project it in there. So r.proj
> does a backward projection for every cell in the new region, and does an
> interpolation to find the most appropriate pixels in the original image
> to re-sample in there.

Maciek / Paul, 

I would add that while gdalwarp uses "forward projection" to establish
the area of the output file in the default case, the actual image
sampling is done via "backward projection" much like r.proj.  

I don't know why the two results are different, but to know whether
it even makes sense to compare you would need to be very careful
in establishing that the two results are for exactly the same region
and resolution.   (Frank reviews the first email, and realizes that
the regions seem to be exact matches with the possible exception
of the odd "Rows: 16746 (16747)" for the output grass region). 

Second, you would want to review very carefully
whether there is anything different about the coordinate systems
that end up getting used.  If one is doing a modest datum shift, and
the other is not, that might well explain differences. 

Thirdly, I would suggest using the -et 0.0 switch with gdalwarp.
Without that, I believe that gdalwarp uses a linear approximation 
for the reprojection as long as the error does not appear to exceed
some fraction of a pixel.  With -et 0.0, gdalwarp will do a full
reprojection operation for each pixel. 

Beyond that, there are various sources of possible mathematical
error that could explain some differences.  However, the differences
in the example you provide do seem to be pretty widespread 
suggestion a substantial issue. 

I would suggest review the "Rows" issue in the GRASS output, 
try again with -et 0.0 and then see where you stand.  If the 
problem persists, then I would suggest reprojecting a grid, 
so that we could careful evaluate some points of difference.   If
you come to the conclusion that gdalwarp is in error, please don't
hesitate to file a bug with the input image, the output image you got
from gdalwarp, and one sample point (pixel) that you believe is in error. 

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent