Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.