[GRASSLIST:6144] r.proj and gdawarp give differernt results

GRASS GRASS user
1

Glynn

I was told that you might know something on this topic and that it was
already disscussed on the grasslist. I have searched the archive but
couldn't find anything relevant. In this situation could you please try to
explain it to me?

The problem is that output of r.proj and gdalwarp is always different,
although the projection parameters and resampling method used are identical
in both approaches.

The difference is that r.proj's product, when compared to gdalwarp's, is
shifted east one resolution unit every few rows.

In case presented here the shift is 14.25 m - resolution of ETM pan. Please
see the attached screendumps. This phenomenon is not limited to the
projections used in this example and occures in any resolution. Depending on
the amount of rotation required during the reprojection, the rows interval,
seperating the identically and differeretly reprojected parts of image,
differs.

Would it be a bug or a result of different reprojection algorithms?

Since these two methods differ which one is right or better?

DETAILS:

    A following ETM pan geotiff:

[pingwin@localhost pingwin]$ gdalinfo p190r025_7p20010524_z33_nn80.tif
Driver: GTiff/GeoTIFF
Size is 17780, 15818
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235629972,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32633"]]
Origin = (501535.875000,5684389.125000)
Pixel Size = (14.25000000,-14.25000000)
Metadata:
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Corner Coordinates:
Upper Left ( 501535.875, 5684389.125) ( 15d 1'19.33"E, 51d18'38.91"N)
Lower Left ( 501535.875, 5458982.625) ( 15d 1'16.03"E, 49d17'0.94"N)
Upper Right ( 754900.875, 5684389.125) ( 18d39'11.29"E, 51d15'13.63"N)
Lower Right ( 754900.875, 5458982.625) ( 18d30'5.38"E, 49d13'49.82"N)
Center ( 628218.375, 5571685.875) ( 16d47'59.05"E, 50d17'0.09"N)
Band 1 Block=17780x1 Type=Byte, ColorInterp=Gray

    ...was gdalwarped into epsg:2180:

[pingwin@localhost pingwin]$ gdalwarp -tr 14.25 14.25 -t_srs epsg:2180
p190r025_7p20010524_z33_nn80.tif
epsg2180_gdalwarp_p190r025_7p20010524_z33_nn80.tif

[pingwin@localhost pingwin]$ gdalinfo
epsg2180_gdalwarp_p190r025_7p20010524_z33_nn80.tif
Driver: GTiff/GeoTIFF
Size is 18600, 16747
Coordinate System is:
PROJCS["ETRS89 / Poland CS92",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",19],
    PARAMETER["scale_factor",0.9993],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",-5300000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2180"]]
Origin = (210751.818018,390197.494682)
Pixel Size = (14.25000000,-14.25000000)
Corner Coordinates:
Upper Left ( 210751.818, 390197.495) ( 14d50'54.92"E, 51d18'17.20"N)
Lower Left ( 210751.818, 151552.745) ( 15d 1'50.83"E, 49d 9'47.62"N)
Upper Right ( 475801.818, 390197.495) ( 18d39'7.96"E, 51d22'40.38"N)
Lower Right ( 475801.818, 151552.745) ( 18d40'3.09"E, 49d13'51.62"N)
Center ( 343276.818, 270875.120) ( 16d47'57.92"E, 50d17'3.16"N)
Band 1 Block=18600x1 Type=Byte, ColorInterp=Gray

    Both geotiffs were then imported into their adequate Grass locations.

    The etm_utm33 location:

GRASS 6.0.cvs:~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name : Universe Transverse Mercator
proj : utm
datum : wgs84
a : 6378137
es : 0.0066943800
zone : 33
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

    The etm_epsg2180 location:

GRASS 6.0.cvs:~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
proj : tmerc
datum : etrs89
a : 6378137
es : 0.0066943800
lat_0 : 0
lon_0 : 19
k : 0.999300
x_0 : 500000
y_0 : -5300000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

    Next the p190r025_7p20010524_z33_nn80 was reprojected from the etm_utm33
    location into etm_epsg2180 location with r.proj:

GRASS 6.0.cvs:~ > r.proj input=p190r025_7p20010524_z33_nn80
location=etm_utm33 mapset=PERMANENT
output=epsg2180_rproj_p190r025_7p20010524_z33_nn80 method=nearest
resolution=14.25

Input Projection Parameters: +proj=utm +no_defs +zone=33 +a=6378137
+rf=298.257223563 +towgs84=0.000,0.000,0.000
Input Unit Factor: 1

Output Projection Parameters: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.999300
+x_0=500000 +y_0=-5300000 +no_defs +a=6378137 +rf=298.257222101
+towgs84=0.000,0.000,0.000
Output Unit Factor: 1
Input:
Cols: 17780 (17780)
Rows: 15818 (15818)
North: 5684389.125000 (5684389.125000)
South: 5458982.625000 (5458982.625000)
West: 501535.875000 (501535.875000)
East: 754900.875000 (754900.875000)
ew-res: 14.250000
ns-res: 14.250000

Output:
Cols: 18600 (18600)
Rows: 16746 (16747)
North: 390197.494682 (390197.494682)
South: 151566.994682 (151552.744682)
West: 210751.818018 (210751.818018)
East: 475801.818018 (475801.818018)
ew-res: 14.250000
ns-res: 14.250000
Allocating memory and reading input map...
Projecting...

After displaying both epsg2180_rproj_p190r025_7p20010524_z33_nn80 and
epsg2180_gdalwarp_p190r025_7p20010524_z33_nn80 it shows that r.proj product
is shifted about one resolution unit (14.25 m in this case) east when
compared with gdalwarp's product. Should both cases be identical?

Maciek

(attachments)

epsg2180_rproj_p190r025_7p20010524_z33_nn80.png
epsg2180_gdalwarp_p190r025_7p20010524_z33_nn80.png

2

Hello Maciek

On Sun, 13 Mar 2005, Maciek Sieczka wrote:

The problem is that output of r.proj and gdalwarp is always different,
although the projection parameters and resampling method used are identical
in both approaches.

[...]

In case presented here the shift is 14.25 m - resolution of ETM pan. Please
see the attached screendumps. This phenomenon is not limited to the
projections used in this example and occures in any resolution. Depending on
the amount of rotation required during the reprojection, the rows interval,
seperating the identically and differeretly reprojected parts of image,
differs.

Would it be a bug or a result of different reprojection algorithms?

Well they are different in the sense that (as far as I understand it) gdalwarp is a forward projection designed for projecting whole images where the image *is* the region of interest and you don't have prior knowledge of the boundary co-ordinates in the new co-ordinate system.

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.

Since these two methods differ which one is right or better?

It depends on whether you are converting a whole image to another projection for unknown future uses or sending it on to someone else (use gdalwarp), or are working on a GIS project covering a certain area and you need to project in data from a different source that covers the region you are working on (use r.proj).

There is plenty of detail in the r.proj man page on how it works but not so much in the gdalwarp page so I could be wrong with what I said about the forward projection. In particular, I'm not too sure at all how it decides what the grid layout in the projected image should be.

But maybe these ideas will put you on the right track to understanding the difference.

Paul

3

From: "Frank Warmerdam" <fwarmerdam@gmail.com>

On Sun, 13 Mar 2005 19:41:30 +0000 (GMT), Paul Kelly
<paul-grass@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).

Sorted out. Was my fault when zooming. See point 6. below that it is
allright now.

Second, you would want to review very carefully
whether there is anything different about the coordinate systems
that end up getting used.

Doublechecked. No differences as to my knowlegde. Both Grass locations
involved where generated from appropriate epsg codes and the same codes are
used with gdalwarp.

If one is doing a modest datum shift, and
the other is not, that might well explain differences.

No datum shift involved since the source projection is utm33/WGS84 and
target is tmerc/GRS80.

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.

Done. The difference if -et is used or not is neglactable and visually
impossible to notice - while difference between r.proj and gdalwarp is very
well visible. However any gdalwarp command I will use from now on will use
the -et 0.0 for the sake of correctness.

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

It persists. Now I will let myself to quote the whole procedure re-done
following your suggestions. All three rasters involved are attached (made
with gdal or exported from grass with r.out.gdal).

    Image "frag_utm33.tif" is my starting point.

    1. reprojecting "frag_utm33.tif" with gdalwarp into
    "frag_gdalwarp_et0_epsg2180.tif":

[pingwin@localhost frag]$ gdalwarp -t_srs epsg:2180 -tr 14.25 14.25 -et 0.0
-co "COMPRESS=DEFLATE" frag_utm33.tif frag_gdalwarp_et0_epsg2180.tif
Creating output file that is 309P x 218L.
:0...10...20...30...40...50...60...70...80...90...100 - done.

[pingwin@localhost frag]$ gdalinfo frag_gdalwarp_et0_epsg2180.tif
Driver: GTiff/GeoTIFF
Size is 309, 218
Coordinate System is:
PROJCS["ETRS89 / Poland CS92",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",19],
    PARAMETER["scale_factor",0.9993],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",-5300000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2180"]]
Origin = (336985.286466,279536.356833)
Pixel Size = (14.25000000,-14.25000000)
Corner Coordinates:
Upper Left ( 336985.286, 279536.357) ( 16d42'26.68"E, 50d21'37.30"N)
Lower Left ( 336985.286, 276429.857) ( 16d42'31.52"E, 50d19'56.77"N)
Upper Right ( 341388.536, 279536.357) ( 16d46'9.43"E, 50d21'41.63"N)
Lower Right ( 341388.536, 276429.857) ( 16d46'14.14"E, 50d20'1.10"N)
Center ( 339186.911, 277983.107) ( 16d44'20.44"E, 50d20'49.22"N)
Band 1 Block=309x26 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)

    2. importing "frag_utm33.tif" into raster "frag_utm33" in the utm33
    Grass location
    3. entering epg:2180 Grass location
    4. importing the "frag_gdalwarp_et0_epsg2180.tif"
    5. reprojecting "frag_utm33" with r.proj:

        a) source utm33 location:

GRASS 6.0.cvs:~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name : Universe Transverse Mercator
proj : utm
datum : wgs84
a : 6378137
es : 0.0066943800
zone : 33
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

        b) target epsg:2180 location

GRASS 6.0.cvs:~ > g.proj -p
-PROJ_INFO-------------------------------------------------
name : Transverse Mercator
proj : tmerc
datum : etrs89
a : 6378137
es : 0.0066943800
lat_0 : 0
lon_0 : 19
k : 0.999300
x_0 : 500000
y_0 : -5300000
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : metre
units : metres
meters : 1

        c) reprojecting

GRASS 6.0.cvs:~ > r.proj input=frag_utm33 location=etm_utm33
mapset=PERMANENT
output=frag_rproj_epsg2180 method=nearest resolution=14.25

Input Projection Parameters: +proj=utm +no_defs +zone=33 +a=6378137
+rf=298.257223563 +towgs84=0.000,0.000,0.000
Input Unit Factor: 1

Output Projection Parameters: +proj=tmerc +lat_0=0 +lon_0=19 +k=0.999300
+x_0=500000 +y_0=-5300000 +no_defs +a=6378137 +rf=298.257222101
+towgs84=0.000,0.000,0.000
Output Unit Factor: 1
Input:
Cols: 299 (299)
Rows: 202 (202)
North: 5580100.500000 (5580100.500000)
South: 5577222.000000 (5577222.000000)
West: 621599.250000 (621599.250000)
East: 625860.000000 (625860.000000)
ew-res: 14.250000
ns-res: 14.250000

Output:
Cols: 309 (720)
Rows: 218 (506)
North: 279542.250000 (281594.250000)
South: 276435.750000 (274383.750000)
West: 336984.000000 (334062.750000)
East: 341387.250000 (344322.750000)
ew-res: 14.250000
ns-res: 14.250000
Allocating memory and reading input map...

Projecting...

    6. both gdalwarp and r.proj product have the *same* number of columns,
    rows and res but the extent and apperance are still somewhat different:

GRASS 6.0.cvs:~ > r.info frag_gdalwarp_et0_epsg2180
+----------------------------------------------------------------------------+
| Layer: frag_gdalwarp_et0_epsg2180 Date: Mon Mar 14 11:55:03 2005|
| Mapset: PERMANENT Login of Creator: pingwin|
| Location: etm_epsg2180| | DataBase: /home/grassdata|
| Title: ( frag_gdalwarp_et0_epsg2180 )|
|----------------------------------------------------------------------------|
||
| Type of Map: raster Number of Categories: 255|
| Data Type: CELL|
| Rows: 218|
| Columns: 309|
| Total Cells: 67362|
| Projection: Transverse Mercator (zone 0)|
| N: 279536.35683273 S: 276429.85683273 Res: 14.25|
| E: 341388.53646595 W: 336985.28646595 Res: 14.25
| Range of data: min = 0 max = 255|
||
| Data Source:|
||
| Data Description:|
| generated by r.in.gdal|
||
+----------------------------------------------------------------------------+

GRASS 6.0.cvs:~ > r.info frag_rproj_epsg2180
+----------------------------------------------------------------------------+
| Layer: frag_rproj_epsg2180 Date: Mon Mar 14 11:44:22 2005|
| Mapset: PERMANENT Login of Creator: pingwin|
| Location: etm_epsg2180|
| DataBase: /home/grassdata|
| Title: ( frag_rproj_epsg2180 )|
|----------------------------------------------------------------------------|
||
| Type of Map: raster Number of Categories: 255|
| Data Type: CELL|
| Rows: 218|
| Columns: 309|
| Total Cells: 67362|
| Projection: Transverse Mercator (zone 0)|
| N: 279542.25 S: 276435.75 Res: 14.25|
| E: 341387.25 W: 336984 Res: 14.25|
| Range of data: min = 14 max = 255|
||
| Data Source:|
||
| Data Description:|
| generated by r.proj|
||
+----------------------------------------------------------------------------+

then I would suggest reprojecting a grid,
so that we could careful evaluate some points of difference.

What resolution, extent and grid interval?

If you come to the conclusion that gdalwarp is in error

What if r.proj is wrong, not gdalwarp? Paul, could that be?

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

Will do.

Maciek

(attachments)

frag_utm33.tif (46 KB)
frag_gdalwarp_et0_epsg2180.tif (48.2 KB)
frag_rproj_epsg2180.tif (73.2 KB)

4

On Mon, 14 Mar 2005 13:39:54 +0100, Maciek Sieczka <werchowyna@epf.pl> wrote:

    6. both gdalwarp and r.proj product have the *same* number of columns,
    rows and res but the extent and apperance are still somewhat different:

Maciek,

I think this is a very significant issue. The difference in
georeferencing seems
to be only a small portion of a pixel, but this can frequently be significant.

To test I tried tracking back one pixel that was different in the two
warpers. I
tracked back the center of a pixel from the GDAL warp output:

Output Location
---------------

Value: 82
Pixel/Line = 39.5, 138.5
Georef = 337548.161466, 277562.731833

Input Location
--------------

cs2cs -f '%.6f' +init=epsg:2180 +to +init=epsg:32633

Georef: 622113.023072 5578151.357107

Pixel/Line: 36.054, 136.782

Value: 82

So, in this case the GDAL result seems to be proper, but the
difference in georeferencing between the r.proj and gdalwarp
output images is approximately 1.28 meters or approximately
0.08882 pixels which would account for the sampling coming
from the neighbouring pixel (which we see in the r.proj result).

I have not traced any other pixels. My suggestion is that you
do what you need to do to get exactly the same georeferencing
in both output products and then compare the results. If there is
still a difference then try tracking back some pixels to see why
they differ (well - at least which is correct).

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

5

Frank

I would be happy to assist you in sorting things out but but based on your
latest reply I don't understand what you want me to do. Please clear it out
to me.

Cheers
Maciek

P.S.
Not being online I usually answer with quite a delay.

From: "Frank Warmerdam" <fwarmerdam@gmail.com>

On Mon, 14 Mar 2005 13:39:54 +0100, Maciek Sieczka <werchowyna@epf.pl>
wrote:

    6. both gdalwarp and r.proj product have the *same* number of
columns,
    rows and res but the extent and apperance are still somewhat
different:

Maciek,

I think this is a very significant issue. The difference in
georeferencing seems
to be only a small portion of a pixel, but this can frequently be
significant.

To test I tried tracking back one pixel that was different in the two
warpers. I
tracked back the center of a pixel from the GDAL warp output:

Output Location
---------------

Value: 82
Pixel/Line = 39.5, 138.5
Georef = 337548.161466, 277562.731833

Input Location
--------------

cs2cs -f '%.6f' +init=epsg:2180 +to +init=epsg:32633

Georef: 622113.023072 5578151.357107

Pixel/Line: 36.054, 136.782

Value: 82

So, in this case the GDAL result seems to be proper, but the
difference in georeferencing between the r.proj and gdalwarp
output images is approximately 1.28 meters or approximately
0.08882 pixels which would account for the sampling coming
from the neighbouring pixel (which we see in the r.proj result).

I have not traced any other pixels. My suggestion is that you
do what you need to do to get exactly the same georeferencing
in both output products and then compare the results. If there is
still a difference then try tracking back some pixels to see why
they differ (well - at least which is correct).

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