3

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
  • Why don't you use -t_srs EPSG:3857. What is the output of gdalinfo if you do that? – pLumo Jun 23 '17 at 10:53
  • $ gdalwarp -t_srs EPSG:3857 LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF out3857.tifCreating output file that is 8902P x 8880L. Processing input file LC08_L1TP_047011_20170518_20170519_01_RT_B4.TIF. 0...10...20...30...40...50...60...70...80...90...100 - done. – dr_times Jun 23 '17 at 10:56
  • 2
    @dr_times Please edit the question, rather than losing important information in comments. – BradHards Jun 23 '17 at 11:08
  • 1
    This is completely solved by using -tr gdalwarp parameter. – xunilk Jun 23 '17 at 11:42
  • 1
    gdalwarp is doing the right thing: keeping the total resolution. Forcing pixelsize with -tr is the wrong way to go for most use-cases. – pLumo Jun 23 '17 at 11:49
10

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

1

That is correct, in Web Mercator the measurement units are not really meters.

See https://stackoverflow.com/a/29672432/4828720

Look at your resulting images, they are probably vastly different in looks.

  • They are in metres. Its projected, and an approximation, but its metres. – BradHards Jun 23 '17 at 11:11
  • 1
    thank you! the document provided on stack overflow is very informative – dr_times Jun 23 '17 at 12:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.