Sunday, 8 March 2015

gdal - Why does reprojection with gdalwarp change pixel size?


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

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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...