13

I have a huge amount of data in the form of CSV containing UTM coordinates as X and Y and an elevation value as Z information. I need to convert these data into a DEM as GeoTIFF for further analysis. In this case, a huge amount means 16 m. lines, with one point in X, Y and Z per line. The points are equally distributed, therefore an interpolation is not needed; each point just needs to be converted into a raster cell.

The original data came without separator, with fixed column widths. I already figured out how to convert the file syntax to use a separator instead of fixed widths and eliminate all space characters, using the stream text editor sed. From here on, normally my workflow would be to import the data into ArcGIS by creating a feature class from the X, Y and Z data and as a second step, converting the point shapefile into a GeoTIFF, using the Point to Raster tool. However, the file I currently have is way too big for this process.

Instead of the above described workflow, I was looking for an efficient alternative and discovered GDAL. However, in gdal_translate, the closest supported format I can find in the supported filetype list, is ASCII grid but no comma separated XYZ. Another difficulty is, that I have UTM coordinates, while most examples seem to use decimal degree coordinates. However, I need to stay within the UTM system (or at least, my output GeoTIFF needs to be in a UTM coordinate system).

So I am looking for a way to convert the CSV XYZ into a GeoTIFF, using GDAL, but so far wasn't able to find examples dealing with this exact problem. I would be very happy for some hints or even code examples.

6
  • Why do you think the GDAL method would be more efficient than the Esri method?
    – artwork21
    Feb 6, 2017 at 13:54
  • The exact example of using a xyz-csv for a tiff is in the documentaion here: gdal.org/gdal_grid.html
    – Matte
    Feb 6, 2017 at 14:06
  • What exactly is the question? Right now the answer is "yes, you can use GDAL to convert". :} Feb 6, 2017 at 14:16
  • The question is how to apply the conversion. Matte's comment seems to provide the solution - I'll try this.
    – Arne
    Feb 6, 2017 at 14:32
  • Ok! Can you provide a minimal example case of data? Do you want an answer in GDAL or would other free tools (eg GMT) be ok as well? Feb 6, 2017 at 15:34

3 Answers 3

19

You can do this using GDAL, it directly supports XYZ format. It doesn't matter if your coordinates are UTM, gdal_translate will output in the same coordinate system.

So to convert to GeoTIFF is as simple as:

gdal_translate test.xyz test.tif

Look at the GeoTIFF doc for output options (such as compression) and the gdal_translate doc for more usage info. In particular, you should specify what the coordinate system is with the -a_srs parameter.

-a_srs srs_def:

Override the projection for the output file. The srs_def may be any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n or a file containing the WKT.

gdal_translate -a_srs EPSG:12345 test.xyz test.tif

Comma/space separated and fixed column widths, with and without a header row are supported.

The supported column separators are space, comma, semicolon and tabulations.

$ head -n 2 test_space.xyz 
x y z
146.360047076550984 -39.0631214488636616 0.627969205379486084

$ gdalinfo test_space.xyz
Driver: XYZ/ASCII Gridded XYZ
Files: test_space.xyz
Size is 84, 66
Coordinate System is `'
Origin = (146.359922066953317,-39.062997159090934)
Pixel Size = (0.000250019195332,-0.000248579545455)
Corner Coordinates:
Upper Left  ( 146.3599221, -39.0629972) 
Lower Left  ( 146.3599221, -39.0794034) 
Upper Right ( 146.3809237, -39.0629972) 
Lower Right ( 146.3809237, -39.0794034) 
Center      ( 146.3704229, -39.0712003) 
Band 1 Block=84x1 Type=Float32, ColorInterp=Undefined
  Min=0.336 Max=0.721 

$ head -n 2 test_commas.xyz 
x, y, z
146.360047076550984, -39.0631214488636616, 0.627969205379486084

$ gdalinfo test_commas.xyz
Driver: XYZ/ASCII Gridded XYZ
etc...

$ head -n 2 test_formatted.xyz
x                       y                       z
146.3600471            -39.06312145             0.627969205

$ gdalinfo test_formatted.xyz
Driver: XYZ/ASCII Gridded XYZ
etc...

The only gotchas I'm aware of are:

  1. The opening of a big dataset can be slow as the driver must scan the whole file to determine the dataset size and spatial resolution; and
  2. The file has to be sorted correctly (by Y, then X).

    Cells with same Y coordinates must be placed on consecutive lines. For a same Y coordinate value, the lines in the dataset must be organized by increasing X values. The value of the Y coordinate can increase or decrease however.

    $ head -n 5 test.csv
    x,y,z
    146.3707979,-39.07778764,0.491866767
    146.3787985,-39.07157315,0.614820838
    146.3637974,-39.07132457,0.555555582
    146.3630473,-39.07579901,0.481217861
    
    $ gdalinfo test.csv
    ERROR 1: Ungridded dataset: At line 3, too many stepY values
    gdalinfo failed - unable to open 'test.csv'.
    
    $ tail -n +2 test.csv| sort -n -t ',' -k2 -k1 > test_sorted.xyz
    
    $ head -n 5 test_sorted.xyz 
    146.3600471,-39.07927912,0.606096148
    146.3602971,-39.07927912,0.603663027
    146.3605471,-39.07927912,0.603663027
    146.3607971,-39.07927912,0.589507282
    146.3610472,-39.07927912,0.581049323
    
    $ gdalinfo test_sorted.xyz
    Driver: XYZ/ASCII Gridded XYZ
    etc...
    
4
  • 3
    I would highly suggest assigning a CRS to the output to make clear what the coordinates are: -a_srs EPSG:12345 Feb 7, 2017 at 10:20
  • 1
    Good point @bugmenot
    – user2856
    Feb 7, 2017 at 10:41
  • where do you run those commands?
    – rocker996
    May 14, 2020 at 14:59
  • In a terminal/command prompt
    – user2856
    May 14, 2020 at 22:30
3

I know this is old, but I just followed the suggestions above and found out gdal_translate is rather slow: 1M lines in xyz format takes dozens of seconds to process - problematic if one needs to process thousands of files containing 1M lines each. Another problem (at least with my files): The .tif's created by gdal_translate had switched coordinates

I have not tested alternatives but just written a Python script using pandas to read the files and rasterio (a gdal wrapper) to write the result. This scripts processes one file in less then 1 second and uses multi-processing for multiple files. A job of 100 files á 1M lines takes 17 seconds on good hardware instead of 50 minutes with a bash loop calling gdal_translate. Until now, the crs is hardcoded to ETRS89/UTM32N (CRS='EPSG:25832'), but can be changed to anything else.

Link to the script: https://gist.github.com/philippkraft/2da0ab4314dd334463fe0e04985bba32

2

Open Data DGM200 of Germany: https://gdz.bkg.bund.de/index.php/default/open-data/digitales-gelandemodell-gitterweite-200-m-dgm200.html

Download the xyz file: https://daten.gdz.bkg.bund.de/produkte/dgm/dgm200/aktuell/dgm200.utm32s.xyzascii.zip

Convert xyz to GeoTIFF file:

gdal_translate -a_srs utm32s.prj dgm200_utm32s.xyz germanyDGM200.tif

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.