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

Maciek Sieczka werchowyna at epf.pl
Mon Mar 14 07:39:54 EST 2005


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

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

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 at 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 at 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 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: frag_utm33.tif
Type: image/tiff
Size: 47097 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20050314/120a199b/frag_utm33.tif
-------------- next part --------------
A non-text attachment was scrubbed...
Name: frag_gdalwarp_et0_epsg2180.tif
Type: image/tiff
Size: 49307 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20050314/120a199b/frag_gdalwarp_et0_epsg2180.tif
-------------- next part --------------
A non-text attachment was scrubbed...
Name: frag_rproj_epsg2180.tif
Type: image/tiff
Size: 74915 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20050314/120a199b/frag_rproj_epsg2180.tif


More information about the grass-user mailing list