[GRASSLIST:6144] r.proj and gdawarp give differernt results
Maciek Sieczka
werchowyna at epf.pl
Sun Mar 13 14:11:10 EST 2005
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 at 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 at 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 at 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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: epsg2180_rproj_p190r025_7p20010524_z33_nn80.png
Type: image/png
Size: 118841 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20050313/515d04db/epsg2180_rproj_p190r025_7p20010524_z33_nn80.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: epsg2180_gdalwarp_p190r025_7p20010524_z33_nn80.png
Type: image/png
Size: 118806 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-user/attachments/20050313/515d04db/epsg2180_gdalwarp_p190r025_7p20010524_z33_nn80.png
More information about the grass-user
mailing list