Sunday, 4 March 2018

gdal - gdal_translate creates images that are mirrored


I'm having a problem with gdal_translate - when translating from an XYZ file to AAIGrid or geotiff the raster is mirrored and "moved" south on one axis.


Mirrored image, lidar data from some open mines


The top part of the attached image is the original file, when displayed in QGis, the bottom file has been translated to AAIGrid format with gdal_translate using


gdal_translate -of AAIGrid -a_srs EPSG:3301 grid.xyz grid.grd

The same happens when I try to translate the XYZ file to Geotiff.


Edit: added file sample. The file looks like this, the first 2 colums are coordinates in EPSG:3301 local coordinate system and the third column is height in metres


719935 6575005 30.709999

719945 6575005 31.08
719955 6575005 30.805
719965 6575005 30.772499
719975 6575005 30.2775
719985 6575005 30.1175
719995 6575005 30.012501
715005 6575015 28.525
715015 6575015 28.715
715025 6575015 28.834999
715035 6575015 28.6875

715045 6575015 28.452499
715055 6575015 29.147499

Output from gdalinfo


gdalinfo 65714_dem_10m.xyz.grid 
Driver: AAIGrid/Arc/Info ASCII Grid
Files: 65714_dem_10m.xyz.grid
65714_dem_10m.xyz.prj
Size is 500, 500
Coordinate System is:

PROJCS["Estonian_Coordinate_System_of_1997",
GEOGCS["GCS_EST97",
DATUM["Estonia_1997",
SPHEROID["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",59.33333333333334],
PARAMETER["standard_parallel_2",58],
PARAMETER["latitude_of_origin",57.51755393055556],

PARAMETER["central_meridian",24],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",6375000],
UNIT["Meter",1]]
Origin = (715000.000000000000000,6575000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Corner Coordinates:
Upper Left ( 715000.000, 6575000.000) ( 27d46'16.00"E, 59d15'32.02"N)
Lower Left ( 715000.000, 6570000.000) ( 27d45'58.29"E, 59d12'50.69"N)
Upper Right ( 720000.000, 6575000.000) ( 27d51'31.03"E, 59d15'22.83"N)

Lower Right ( 720000.000, 6570000.000) ( 27d51'12.91"E, 59d12'41.51"N)
Center ( 717500.000, 6572500.000) ( 27d48'44.56"E, 59d14' 6.79"N)
Band 1 Block=500x1 Type=Float32, ColorInterp=Undefined

Answer



The .xyz file is sorted in an unusual way, X ascending and Y ascending. Most GDAL drivers can handle that orientation by flipping the image internally. Thus the "upper" and "lower" coordinates are flipped too:


Driver: XYZ/ASCII Gridded XYZ
Files: test.xyz
Size is 500, 2
Coordinate System is `'
Origin = (715000.000000000000000,6575000.000000000000000)

Pixel Size = (10.000000000000000,10.000000000000000)
Corner Coordinates:
Upper Left ( 715000.000, 6575000.000)
Lower Left ( 715000.000, 6575020.000)
Upper Right ( 720000.000, 6575000.000)
Lower Right ( 720000.000, 6575020.000)
Center ( 717500.000, 6575010.000)
Band 1 Block=500x1 Type=Float32, ColorInterp=Undefined
Min=28.452 Max=31.080
NoData Value=0


If you sort the .xyz file by X ascending and Y descending, the metadata appears in a logical way:


Driver: XYZ/ASCII Gridded XYZ
Files: testsort.xyz
Size is 500, 2
Coordinate System is `'
Origin = (715000.000000000000000,6575020.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Corner Coordinates:
Upper Left ( 715000.000, 6575020.000)

Lower Left ( 715000.000, 6575000.000)
Upper Right ( 720000.000, 6575020.000)
Lower Right ( 720000.000, 6575000.000)
Center ( 717500.000, 6575010.000)
Band 1 Block=500x1 Type=Float32, ColorInterp=Undefined
Min=28.452 Max=31.080
NoData Value=0

Unfortunately, the AAIGrid format is not designed for the first orientation, and noone has implemented to catch this issue by an error message or internally flipping back the coordinates. As a consequence, the header for AAIGrid gets generated wrong:


ncols        500

nrows 2
xllcorner 715000.000000000000
yllcorner 6574980.000000000000
cellsize 10.000000000000
NODATA_value 0

The following data is always linewise from upper left to lower right.


As a kind of hack, you can use gdalwarp instead of gdal_translate to save your raster as tif, and then translate to AAIGrid:


gdalwarp test.xyz 3301.tif
gdalinfo 3301.tif

gdal_translate -of AAIGRID 3301.tif 3301.asc

which delivers the correct header:


ncols        500
nrows 2
xllcorner 715000.000000000000
yllcorner 6575000.000000000000
cellsize 10.000000000000
NODATA_value 0

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...