OSGeo
gdal_merge.py doesn't support bottom-to-top images
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
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.
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!
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`
