GDAL – Why Does Reprojection with gdalwarp Change Pixel Size?

coordinate systemgdalgdalwarp

Is it okay that reprojection from UTM (EPSG:32612) to Web Mercator (EPSG:3857) with gdalwarp changes pixel size? Why does this happen? Does that mean that unit is not meters? If so, than how to convert pixel size from UTM units to Web Mercator units?

Command:

gdalwarp -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs" LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF output.tif

Input image is Landsat 8 from AWS:
http://landsat-pds.s3.amazonaws.com/c1/L8/047/011/LC08_L1TP_047011_20170518_20170519_01_RT/LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF

Original image info:

$ gdalinfo LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF
Driver: GTiff/GeoTIFF
Files: /home/user/Desktop/Landsat_sample/LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF
Size is 8451, 8491
Coordinate System is:
PROJCS["WGS 84 / UTM zone 12N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-111],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32612"]]
Origin = (488985.000000000000000,7851615.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  488985.000, 7851615.000) (111d17'58.67"W, 70d46' 6.92"N)
Lower Left  (  488985.000, 7596885.000) (111d16' 8.84"W, 68d29' 3.48"N)
Upper Right (  742515.000, 7851615.000) (104d25'49.46"W, 70d39' 3.78"N)
Lower Right (  742515.000, 7596885.000) (105d 5'39.63"W, 68d22'48.67"N)
Center      (  615750.000, 7724250.000) (108d 1'24.05"W, 69d36' 5.29"N)
Band 1 Block=512x512 Type=UInt16, ColorInterp=Gray

Reprojected image info:

$ gdalinfo output.tif
Driver: GTiff/GeoTIFF
Files: out.tif
       out.tif.aux.xml
Size is 8902, 8880
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-12389818.176861140877008,11323689.164804801344872)
Pixel Size = (85.895238477817642,-85.895238477817642)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-12389818.177,11323689.165) (111d17'58.67"W, 70d46' 7.79"N)
Lower Left  (-12389818.177,10560939.447) (111d17'58.67"W, 68d22'48.98"N)
Upper Right (-11625178.764,11323689.165) (104d25'50.73"W, 70d46' 7.79"N)
Lower Right (-11625178.764,10560939.447) (104d25'50.73"W, 68d22'48.98"N)
Center      (-12007498.470,10942314.306) (107d51'54.70"W, 69d36'28.82"N)
Band 1 Block=8902x1 Type=UInt16, ColorInterp=Gray
  Min=0.000 Max=52607.000 
  Minimum=0.000, Maximum=52607.000, Mean=17537.063, StdDev=16780.853
  Metadata:
    STATISTICS_MAXIMUM=52607
    STATISTICS_MEAN=17537.062597405
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=16780.852772678

Best Answer

gdalwarp is doing the right thing: preserving total resolution of your image by changing the pixel-size.

WGS 84 / Pseudo-Mercator projection is heavily distorted when moving away from the equator. Thus, it could be discussed if the units should be called "Pseudo-meters". One meter in reality is approximately 1/cos(lat) pseudo-meters.

You can calculate the approximate pixel size for your raster image:

30m / cos(69.53°) = 85.78 pseudo-meters

Compare this with your gdalinfo output.

Tissot Indicatrix

CC BY-SA 3.0, Author: Stefan Kühn

Related Question