Bilinear Interpolation Incorrect

[johanstenberg]

Hi,

My name is Johan Stenberg and I work as a software consultant for GeoTrellis. When looking for a reference implementation of bilinear interpolation, I looked at GDAL's bilinear interpolation and I think it's broken.

I have here attached a file for which the command :

"gdalwarp -r bilinear -s_srs '+proj=longlat +datum=WGS84 +no_defs' -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' nlcd_tile_wsg84.tif nlcd_tile_webmercator-bilinear.tif"

Generates a bad file.

It's bad because in the resulting file (there might be more errors than this one) one value is -31488.0. There is no value in the original file that is lower than -30000, so that should not be able to happen in regular bilinear interpolation.

I've been trying to spot the error myself but I haven't had any success.

Thanks,

Johan

[johanstenberg]

TIFF file used in bug report.

[Even Rouault]

This is mathematically impossible. The data type of the file is UInt16, so negative values cannot happen. Below the statistics on the source and target files (with trunk, very slightly different with 1.11 due to small changes in trunk). The min/max values are the same.

$ gdalinfo -stats nlcd_tile_wsg84.tif
Driver: GTiff/GeoTIFF
Files: nlcd_tile_wsg84.tif
Size is 512, 512
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-96.952496999999994,38.143939000000003)
Pixel Size = (0.000325872016369,-0.000325872016369)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -96.9524970,  38.1439390) ( 96d57' 8.99"W, 38d 8'38.18"N)
Lower Left  ( -96.9524970,  37.9770925) ( 96d57' 8.99"W, 37d58'37.53"N)
Upper Right ( -96.7856505,  38.1439390) ( 96d47' 8.34"W, 38d 8'38.18"N)
Lower Right ( -96.7856505,  37.9770925) ( 96d47' 8.34"W, 37d58'37.53"N)
Center      ( -96.8690738,  38.0605158) ( 96d52' 8.67"W, 38d 3'37.86"N)
Band 1 Block=512x8 Type=UInt16, ColorInterp=Gray
  Minimum=71.000, Maximum=237.000, Mean=210.668, StdDev=28.888
  Metadata:
    STATISTICS_MAXIMUM=237
    STATISTICS_MEAN=210.66777801514
    STATISTICS_MINIMUM=71
    STATISTICS_STDDEV=28.888394822497


$ gdalinfo -stats nlcd_tile_webmercator-bilinear.tif
Driver: GTiff/GeoTIFF
Files: nlcd_tile_webmercator-bilinear.tif
Size is 448, 569
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 = (-10792702.597176384180784,4599779.568733124993742)
Pixel Size = (41.464747379624953,-41.464747379624953)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-10792702.597, 4599779.569) ( 96d57' 8.99"W, 38d 8'38.18"N)
Lower Left  (-10792702.597, 4576186.127) ( 96d57' 8.99"W, 37d58'37.43"N)
Upper Right (-10774126.390, 4599779.569) ( 96d47' 8.25"W, 38d 8'38.18"N)
Lower Right (-10774126.390, 4576186.127) ( 96d47' 8.25"W, 37d58'37.43"N)
Center      (-10783414.494, 4587982.848) ( 96d52' 8.62"W, 38d 3'37.98"N)
Band 1 Block=448x9 Type=UInt16, ColorInterp=Gray
  Minimum=71.000, Maximum=237.000, Mean=210.680, StdDev=26.500
  Metadata:
    STATISTICS_MAXIMUM=237
    STATISTICS_MEAN=210.67953254456
    STATISTICS_MINIMUM=71
    STATISTICS_STDDEV=26.499834325815