Opened 8 years ago

Closed 8 years ago

#6320 closed defect (invalid)

gdalwarp output resolution inconsistant with the one of other software packages

Reported by: qunidaye Owned by: warmerdam
Priority: normal Milestone:
Component: default Version: 2.0.0
Severity: normal Keywords: Sinusoidal projection
Cc:

Description

When I use gdalwarp.exe of GDAL 2.0 to convert an raster with Sinusodial projection to latlong, the resolution of the result data always increase much. Below are the test example and data.

Command: gdalwarp.exe -t_srs "+proj=latlong" in.img out.img

the gdalinfo.exe result of the in.img is

Size is 2400, 2400
Coordinate System is:
PROJCS["Sinusoidal",
    GEOGCS["GCS_Unknown",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["longitude_of_center",0.0],
    UNIT["Meter",1]]
Origin = (4447802.078700000400000,4447802.078700000400000)
Pixel Size = (463.312716529999990,-463.312716529999990)
Metadata:
  Band_1=band1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 4447802.079, 4447802.079) ( 52d12'40.66"E, 40d 9'52.45"N)
Lower Left  ( 4447802.079, 3335851.559) ( 46d 9'49.03"E, 30d 8'31.10"N)
Upper Right ( 5559752.598, 4447802.079) ( 65d15'50.82"E, 40d 9'52.45"N)
Lower Right ( 5559752.598, 3335851.559) ( 57d42'16.28"E, 30d 8'31.10"N)
Center      ( 5003777.339, 3891826.819) ( 54d55' 1.11"E, 35d 9'19.23"N)
Band 1 Block=2400x1 Type=Byte, ColorInterp=Undefined
  Description = band1

the gdalinfo.exe result of the out.img is

Size is 5672, 2976
Coordinate System is:
GEOGCS["unnamed ellipse",
    DATUM["unknown",
        SPHEROID["unnamed",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (46.163618708480982,40.164568365792846)
Pixel Size = (0.003367378242600,-0.003367378242600)
Metadata:
  AREA_OR_POINT=Area
  Band_1=band1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  46.1636187,  40.1645684) ( 46d 9'49.03"E, 40d 9'52.45"N)
Lower Left  (  46.1636187,  30.1432507) ( 46d 9'49.03"E, 30d 8'35.70"N)
Upper Right (  65.2633881,  40.1645684) ( 65d15'48.20"E, 40d 9'52.45"N)
Lower Right (  65.2633881,  30.1432507) ( 65d15'48.20"E, 30d 8'35.70"N)
Center      (  55.7135034,  35.1539095) ( 55d42'48.61"E, 35d 9'14.07"N)
Band 1 Block=5672x1 Type=Byte, ColorInt
erp=Gray
  Description = band1

As can be see, the final resolution of the pixel is about 0.003367378242600, but in ArcGIS and ENVI is about 0.005.

By the way, I have try to use the c++ to re-project the projection, but it also with the same result (I have test it both in version of 1.11 and 2.0).

Attachments (2)

in.rar (194.0 KB ) - added by qunidaye 8 years ago.
in.img
out.rar (284.1 KB ) - added by qunidaye 8 years ago.
out.tif

Download all attachments as: .zip

Change History (5)

by qunidaye, 8 years ago

Attachment: in.rar added

in.img

by qunidaye, 8 years ago

Attachment: out.rar added

out.tif

comment:1 by Even Rouault, 8 years ago

Summary: problem of convert Sinusoidal projection to other projectionsgdalwarp output resolution inconsistant with the one of other software packages

The logics to guess the output resolution must be different between packages. You can always force the target resolution with the -tr option.

The logic used by GDAL is the following one :

  • transform the coordinates of the upper-left and lower-right corners to the target projection
  • divide this distance by the number of pixels of the diagonal of the input image

Demo :

$ gdaltransform -s_srs "+proj=sinu" -t_srs "+proj=longlat"
4447802.079 4447802.079
52.2112937270165 40.1645683684946 0

5559752.598 3335851.559
57.7045233815621 30.1419724337471 0

$ python
>>> import math
>>> math.sqrt((57.7045233815621-52.2112937270165)**2 + (30.1419724337471-40.1645683684946)**2)
11.429260750757077
>>> math.sqrt(2400**2 + 2400 ** 2)
3394.1125496954282
>>> 11.429260750757077 / 3394.1125496954282
0.0033673782420039918

When tring with the other diagonal (lower-left and upper-right), this would give:

$ gdaltransform -s_srs "+proj=sinu" -t_srs "+proj=longlat"
4447802.079, 3335851.559  
46.1636187114771 30.1419724337471 0

5559752.598, 4447802.079
65.2641171499666 40.1645683684946 0

$ python
>>> math.sqrt((46.1636187114771-65.2641171499666)**2 + (30.1419724337471-40.1645683684946)**2)
21.570384091850492
>>> 21.570384091850492 / 3394.1125496954282
0.0063552353600607962

And if you average 0.0033673782420039918 + 0.0063552353600607962) / 2, you get 0.0048 . Perhaps this is the logic they use.

comment:2 by qunidaye, 8 years ago

Thank you.

comment:3 by Jukka Rahkonen, 8 years ago

Resolution: invalid
Status: newclosed
Note: See TracTickets for help on using tickets.