gdal_merge.py doesn't support bottom-to-top images

pcace commented 9 months ago

Expected behavior and actual behavior.

i convert XYZ files to GeoTIFF:

gdal_translate ${xyzFile} ${xyzFile}.tif -a_srs ${srs} and then merge them: gdal_merge.py -o merged.tif *.tif the result is always:

gdal_merge.py -o merged.tif *.tif
ERROR 1: Attempt to create 3000x-999 dataset is illegal,sizes must be larger than zero.
Creation failed, terminating gdal_merge.

expected would be a proper merged geoTIFF. this only happens if i use the GeoTIFFS created by gdal_translate. using other GeoTIFFs the merge command works

Steps to reproduce the problem.

Here are the used XYZ Files. they are Projected in EPSG:25832. so we use: gdal_translate dgm1_32_568_5941_1_hh.xyz dgm1_32_568_5941_1_hh.xyz.tif -a_srs EPSG:25832`

dgm1_32_568_5941_1_hh.xyz.gz dgm1_32_570_5941_1_hh.xyz.gz dgm1_32_570_5942_1_hh.xyz.gz

and to merge gdal_merge.py -o merged.tif *.tif

Operating system

POPOS22.04 & Debian 12

GDAL version and provenance

Version: GDAL 3.4.3, released 2022/04/22

GDALINFO

the gdalinfo on on resulting tiff file to merge looks like this:

Driver: GTiff/GeoTIFF
Files: dgm1_32_570_5941_1_hh.xyz.tif
Size is 1000, 1000
Coordinate System is:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
        BBOX[38.76,6,84.33,12]],
    ID["EPSG",25832]]
Data axis to CRS axis mapping: 1,2
Origin = (569999.500000000000000,5940999.500000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  569999.500, 5940999.500) ( 10d 3'29.23"E, 53d36'47.25"N)
Lower Left  (  569999.500, 5941999.500) ( 10d 3'30.04"E, 53d37'19.60"N)
Upper Right (  570999.500, 5940999.500) ( 10d 4'23.63"E, 53d36'46.76"N)
Lower Right (  570999.500, 5941999.500) ( 10d 4'24.45"E, 53d37'19.12"N)
Center      (  570499.500, 5941499.500) ( 10d 3'56.84"E, 53d37' 3.18"N)
Band 1 Block=1000x2 Type=Float32, ColorInterp=Gray

jratike80 commented 9 months ago

The problem comes probably from the interpretation of the pixel size and the orientation of the image. Have a look at the source data

gdalinfo dgm1_32_568_5941_1_hh.xyz
...
Origin = (567999.500000000000000,5940999.500000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Corner Coordinates:
Upper Left  (  567999.500, 5940999.500)
Lower Left  (  567999.500, 5941999.500)
Upper Right (  568999.500, 5940999.500)
Lower Right (  568999.500, 5941999.500)
...

How to read this information:

the origin is at the top-left corner
the pixel size is 1,1 so each step to the right is increasing X and each step towards bottom is increasing Y
the top Y is 5940999.500 and the bottom Y is 5941999.500 that matches with the positive Y pixel size
but in the EPSG:25832 system coordinate Y does not advance from top to bottom but from bottom to top, from the equator towards the north pole

That means that for showing the tiff image in EPSG:25832 the image must be flipped. See how the TIFF image looks with IrfanView that show the raw image, and in QGIS that knows to flip it.

Maybe gdal_merge.py cannot handle flipped images. Gdalbuildvrt cannot either

gdalbuildvrt merge.vrt *.tif
0...10...20...30.Warning 1: gdalbuildvrt does not support positive NS resolution. Skipping dgm1_32_568_5941_1_hh.tif

I did not make a proper test but I believe that you will have success by using gdalwarp that turns the image into north-up

gdalwarp -t_srs epsg:25832 dgm1_32_568_5941_1_hh.xyz dgm1_32_568_5941_1_hh.tif

I cannot say if there is a bug and where it actually could be, but I do believe that is it not easy for an average GDAL user to find out what happens. At least gdalbuildvrt gives a more meaningful error message.

pcace commented 9 months ago

Hi, thank you so much for your help! I think gdalwarp does the trick for me. The question is, if this is now a bug and should be fixed somehow, or can this be closed as this is kind of expected behaviour? Please close it here if you mean this is not a bug. i am not sure...

Thanks anyways for your extremely fast help!

becryan commented 5 months ago

As pointed out as well above by @jratike80 , -- 'the pixel size is 1,1 so each step to the right is increasing X and each step towards bottom is increasing Y'

You can also remedy the issue upstream in your gdal_translate or gdalwarp using the '-tr' parameter which defines pixel size -- as '1 -1' so that you have increasing x, decreasing y. e.g:

gdal_translate ${xyzFile} ${xyzFile}.tif -a_srs ${srs} -tr 1.0 -1.0`

OSGeo / gdal