Pixel offsets occur using gdal_merge.py with -ul_lr argument

[ken.boss@…]

If I run gdal_merge.py with a -ul_lr argument, output pixels are shifted left one pixel and and up one pixel from where they were in the input image. Without a -ul_lr argument, there is no shift in the output.

A simple sample command line that produces shifted output is:

gdal_merge.py -o out.img -v -f HFA -ul_lr 589177 5291586 590519 5290159 in.img

[ken.boss@…]

After further investigation, I need to recharacterize the problem. I am using
the ERDAS HFA file format, and suspect my problem lies there (more below). 
Image UL LR coordinates given here are as reported by ERDAS Imagine.

Running gdal_merge.py without a -ul_lr argument does not produce shift in output
pixels, as stated earlier.  For example,

gdal_merge.py -o no_ullr.img -v -f HFA test.img

gives me back exactly what I put in (ignoring header info); that is, it produces
a copy of my input image.

Running gdal_merge.py with a -ul_lr argument produces shifts in output pixels,
but in a different manner than I described previously.  I'll use what amounts to
a subsetting operation for simplicity's sake, e.g.,

gdal_merge.py -o straight_ullr.img -v -f HFA -ul_lr 590475 5291461 590683
5291255 test.img

produces an output image with these corner coordinates:

ulx 590475.5
uly 5291460.5
lrx 590682.5
lry 5291255.5

Essentially, the output image has contracted by half a unit (pixel, meter, same
thing in this case) in all directions relative to what I requested. "Missing
pixels" were taken from the center of the output image, both horizontally and
vertically, with the effect that pixels in the each quadrant of the output image
have shifted half a pixel toward the center of the image.

In order to get what I actually wanted out of the script, I can fudge the
coordinates of my -ul_lr argument, like:

gdal_merge.py -o tweak_ullr.img -v -f HFA -ul_lr 590474.5 5291461.5 590683.5
5291254.5 test.img

which produces an output image with no missing or shifted pixels and these
corner coordinates:

ulx 590475 
uly 5291461 
lrx 590683 
lry 5291255

I suspect that all of this relates to the fact that ERDAS takes the center of a
pixel to be it's coordinate location, while GDAL appears to want to take a
corner of a pixel to be the coordinate location, and something is getting lost
in the translation...

[ken.boss@…]

Converted my test.img input to geotiff (test.tif) and ran:

gdal_merge.py -o straight_ullr.tif -v -f GTiff -ul_lr 590475 5291461 590683
5291255 test.tif

Results were the same as when I used the HFA format, so my suspicions were
incorrect and this bug does not in fact appear to relate to the file format.

[warmerdam]

This should be reviewed before the 1.2.0 release.

[warmerdam]

Ken,

I realize it has been a long time since you filed this.   Sorry for the delay
in digging in.  

I used this file as input:
warmerda@gdal2200[150]% gdalinfo utm.tif 
Driver: GTiff/GeoTIFF
Size is 512, 512
Coordinate System is:
PROJCS["NAD27 / UTM zone 11N",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.9786982138982,
                AUTHORITY["EPSG","7008"]],
            AUTHORITY["EPSG","6267"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4267"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26711"]]
Origin = (440720.000000,3751320.000000)
Pixel Size = (60.00000000,-60.00000000)
Corner Coordinates:
Upper Left  (  440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N)
Lower Left  (  440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)
Upper Right (  471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N)
Lower Right (  471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)
Center      (  456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N)
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray

I ran this command:
gdal_merge.py -o out.tif -v utm.tif -ul_lr 440780 3751260 443780 3748260

I got this output:
warmerda@gdal2200[149]% gdalinfo out.tif 
Driver: GTiff/GeoTIFF
Size is 50, 50
Coordinate System is `'
Origin = (440780.000000,3751260.000000)
Pixel Size = (60.00000000,-60.00000000)
Corner Coordinates:
Upper Left  (  440780.000, 3751260.000) 
Lower Left  (  440780.000, 3748260.000) 
Upper Right (  443780.000, 3751260.000) 
Lower Right (  443780.000, 3748260.000) 
Center      (  442280.000, 3749760.000) 
Band 1 Block=50x50 Type=Byte, ColorInterp=Gray

I also overlaid the two images in a viewer and everything matched up to
a subpixel level. 

I suspect your request was made on non-pixel boundaries.  I believe
gdal_merge.py computes a window in the source and destination and then
computes the output bounds based on that window rather than based on the
request made.  

The pixels boundaries likely are on half meters amounts due to the middle/edge
pixel differences between erdas format and the GDAL data model. 

If you don't have any more information suggesting a bug I will close this
in a week or two. 

[neteler@…]

Frank,

> The pixels boundaries likely are on half meters amounts due to the middle/edge
> pixel differences between erdas format and the GDAL data model.

Is this an exeption or could this be somehow treated?

From the GDAL documentation I didn't found out yet if cell centers or true pixel
corners are displayed in gdalinfo, but maybe my oversight.

Best regards,

 Markus

[warmerdam]

Markus, 

It is inevitable that gdal_merge.py needs to map a pixel window to a
pixel window.  So, it is just being honest about what it really does when
it sets output bounds slightly different from what was requested.

For simple subsetting, I think it is much better to use gdalwarp which 
can do interpolation, and doesn't have alot of the built-in limitations that
gdal_merge.py does.  Gdal_merge.py was really intended to be a quick hack to
demonstrate the capability of the python interface. 

[ken.boss@…]

Frank--

Granted gdalwarp would be better for simple subsets.  I used a subset operation 
as an example to keep it simple; in fact I am using gdal_merge.py for automated 
production mosaicking (thanks!).

Anyway, since I've come to understand the problem (pixel corners vs. pixel 
centers), I'd no longer consider this a bug.  As Markus points out, though, it 
would be useful if the documentation (perhaps the documentation of the 
utilities' command-line args) were to state explicitly that GDAL takes bounding 
box coordinates to be the outer corners of corner pixels.  I don't expect that 
there's any means of addressing this programmatically; users simply need to be 
aware of what bounding box language they're speaking.

[neteler@…]

Frank,

please let me express my high interest in understanding the
problem of pixel center registration and pixel boundary
coordinates. To be precise, it would be great to have a
FAQ entry (or whereever), which contains something like this:

%<---------------- suggestion -------------

</p><p>
</p></li><li> <b>GDAL and pixel coordinates</b>
<p>

</p><p>
<ul>
<li> gdalinfo reports 'Corner Coordinates' which is the boundary
  box of the raster matrix</li>
<li> gdalinfo reports 'Center' coordinates as center of the central
  pixel of the raster matrix</li>
<li> gdal_translate/gdalwarp use center/corner coordinates depending on 
  the input/output format specifications.</li>
</ul>
</p><p>
</p>

%<---------------- end suggestion (please correct if wrong) -------------

Currently I have a hard time to convince a user that my script [1] to
create GeoTIFF from SRTM BIL files works correctly. IMHO yes, but
he is confused about how GDAL works.

[1] http://mpa.itc.it/rs/srtm/index.html , at page bottom

We would really appreciate the clarification of this problem.

Thanks for your time,

 Markus Neteler

[Jukka Rahkonen]

Issue was created only because of the confuse from one software reporting the image extents measured from centers of corner pixels, another measures from the out-most corners of corner pixels.