Bug in gdalwarp's simple nearest neighbor? Or more likely I need some help

Question

I'm trying to understand how to use the GDAL utilities a little better. I think I've defined something incorrectly or I've found a bug in gdalwarp. My bet is that I've defined something incorrectly, I would love some help in understanding this better. I'm getting answers for a simple nearest neighbor reprojection that I don't understand. I apologize in advance for the length of this post, but it's hard to describe in any fewer words. If someone who's a gdal master could follow it along and show me where I fell of the path, I'd really appreciate it.

I start with a 5x5 raw 2-byte image, with values 1 through 25 ( [[1, 2, 3, 4, 5][6,...,10][...]]) file:ease_subset_324_376_5x5.bin

I'm letting this represent specific grid cells in a 25km EASE2 projection.
ISPRS Int. J. Geo-Inf. 2012, 1, 32-45; doi:10.3390/ijgi1010032
This projection has a proj4 string:

'+proj=laea +lat_0=90. +lon_0=0. +x_0=0. +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m'

In my case the upper left corner of the upper left grid cell is located at:
-900000.00, -400000.00

Which means the lower right corner of the lower right grid cell is located at:
-900000.00 + 5 * 25000., -400000.00 - (5 * 25000.) or -775000., -525000.

I want to convert this to a format that contains the projection and location information so I create a VRT file to desibe the raw data:

<VRTDataset rasterXSize="5" rasterYSize="5">
  <VRTRasterBand dataType="Int16" band="1" subClass="VRTRawRasterBand">
    <SourceFilename relativetoVRT="1">ease_subset_324_376_5x5.bin</SourceFilename>
    <ImageOffset>0</ImageOffset>
    <PixelOffset>2</PixelOffset>
    <LineOffset>10</LineOffset>
    <ByteOrder>LSB</ByteOrder>
  </VRTRasterBand>
</VRTDataset>

And use gdal_translate to create a GeoTiff that describes my data:

gdal_translate subset.vrt subset_EASE2_N25km_gdal.tif \
    -a_srs '+proj=laea +lat_0=90. +lon_0=0. +x_0=0. +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m' \
    -a_ullr -900000.00 -400000.00 -775000.00 -525000.00

Next I want to reproject this image to a 114x105 2km polar stereographic grid. The upper left corner of the upper left gridcell is located: -366000.00 -782000.00 and the lower right corner of the lower right gridcell is: -366000 + (114 * 2000.), -782000.00 - (105 * 2000.) or: -138000. -992000.

I use gdal warp to project my image into my new stereogrphic projection.

gdalwarp -t_srs "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" \
     -te -366000.00 -992000.00  -138000.00 -782000.00 \
     -tr 2000 2000 subset_EASE2_N25km_gdal.tif subset_Moge2km_gdal.tif

I now have an image 114x105 with values 0-25 like I expect, but there are a few places, where I don't understand the values in the reprojected data.

For example: The value in gridcell[50, 26] => 5.

I would presume nearest neighbor sampling, would take the center of that gridcell, project it into the original image and select the closest value from there, but that's not what I see.

the location of the center of polarstereo gridcell[50, 26] is (2km cells)
-366000.00 + 2000./2 + 50 * 2000., -782000.00 - 2000./2 - 26 * 2000.
or: -265000.00, -835000.00

Using gdaltransform I can find the location in my original image:

 gdaltransform
      -s_srs "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
      -t_srs "+proj=laea +lat_0=90. +lon_0=0. +x_0=0. +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m"

-265000.00, -835000.00
-800001.11880625 -414546.034290511 0

Converting this x/y onto my original EASEgrid2 grid cell.

Center of upper left gridcell = -900000 + 25000./2, -400000 - 25000./2. or: -887500., -412500.

xmeters = x - x0 = -800001.11880625 - (-887500.) = 87498.88119374996
ymeters = y0 - y = -412500 - (-414546.034290511) =  2046.034290510986

Convert to gridcells by dividing by resolution:

x_gridcell = 87498.88119374996 / 25000. =  3.4999552
y_gridcell = 2046.034290510986 / 25000. =  0.081841372

This to me would say that I would want the value at 3,0, but I was getting the value at 4,0. grid[3,0] = 4, grid[4,0] = 5

Also if you just look at the meters:

back_projected_point = -800001.11880625, -414546.034290511
grid at 3,0 = g3 => -812500.00      -412500.00
grid at 4,0 = g4 => -787500.00      -412500.00

distance(back_projected_point, g3) = 12665.239
distance(back_projected_point, g4) = 12667.447

Again, this to me looks like I should be getting the value at gridcell[3,0], but that's not what I see.

Following this logic: All of these points in my polar stereo projected image are incorrect.

points =  [ [ 47, 88 ], $
            [ 59, 87 ], $
            [ 60, 88 ], $
            [ 61, 89 ], $
            [ 62, 90 ], $
            [ 50, 26 ], $
            [ 49, 25 ], $
            [ 48, 24 ], $
            [ 47, 23 ], $
            [ 46, 22 ], $
            [ 45, 21 ] ]

Again, all I can think of is that I'm not properly specifying the gdal parameters, or there's something wrong with gdal.

Answer 1

Apparently, I am not specifying the gdal parameters correctly.

After getting bitten by this almost 2 years later in a different manner, kwbeam answered a question regarding horizontal artifacts in gdal which applies directly to this question as well.

The short answer is that you have to make the error tolerance very small in the gdalwarp command. I presume there are some optimizations that cause these reprojections to provide unexpected answers in gdalwarp, but by specifying the -et {tol} to be very small you can get the expected answer.

so by changing my warp command to this:

gdalwarp -t_srs "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 \ 
 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" \
 -te -366000.00 -992000.00  -138000.00 -782000.00 \
 -et .001 \
 -tr 2000 2000 subset_EASE2_N25km_gdal.tif subset_Moge2km_gdal.tif

I now get the output I expected for all the 11 pixels above.