Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have just exported a tif DEM to XYZ format with SAGA GIS. When I check this with gdalinfo it states that the format is unsupported?

Seeing as I converted this with saga gis which used gdal I find this rather strange.

Can anyone suggest what I did wrong?

thanks.

share|improve this question
add comment

2 Answers

Since you exported a DEM, I'm assuming you had gridded data in your XYZ, and it's reasonable to expect that gdal_translate would have handled it. I'm guessing that the problem is that SAGA exports XYZ with only the good points from your grid, and not the missing data points. It seems that gdal_translate needs all the points from the grid in order to work. I wrote this little python script to convert some similar XYZ data. I haven't made this very flexible, but I'm sure you could modify it to fit your needs. The Ipython Notebook version is here: http://nbviewer.ipython.org/urls/raw.github.com/rsignell-usgs/ipython-notebooks/master/xyz2tif.ipynb, and the plain ol' python code is below:

#!/usr/bin/env python 
# xyz2tif.py
# convert gridded XYZ data (but without missing values) to GeoTIFF
# data looks like this:
"""
379528.00 4577221.00 -5.38
379529.00 4577221.00 -5.38
379530.00 4577221.00 -5.40
379531.00 4577221.00 -5.41
379515.00 4577222.00 -5.34
379516.00 4577222.00 -5.36
"""

from pylab import *
from osgeo import osr,gdal

# UTM ZONE 19, NAD83 data
data =genfromtxt('/rps/bathy/muskegat/EW_Survey1.xyz')
x=data[:,0]
y=data[:,1]
z=data[:,2]

# find the extent and difference 
xmin=min(x);xmax=max(x)
ymin=min(y);ymax=max(y)
mdx=abs(diff(x))
mdy=abs(diff(y))

# determine dx and dy from the median of all the non-zero difference values
dx=median(mdx[where(mdx>0.0)[0]])
dy=median(mdy[where(mdy>0.0)[0]])

#construct x,y,z of complete grid
xi=arange(xmin,xmax+dx,dx)
yi=arange(ymin,ymax+dy,dy)
zi=ones((len(yi),len(xi)))*NaN
shape(zi)

# calculate indices in full grid (zi) to stick the input z values
ix=numpy.round((x-xmin)/dx).astype(int)
iy=numpy.round((y-ymin)/dy).astype(int)
zi[iy,ix]=z

zi=flipud(zi)

# write as 32-bit GeoTIFF using GDAL
ny,nx = zi.shape
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create('output.tif', nx, ny, 1, gdal.GDT_Float32)

# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
ds.SetGeoTransform( [ xmin, dx, 0, ymax, 0, dy ] )

# set the reference info 
srs = osr.SpatialReference()

# UTM zone 19, North=1
srs.SetUTM(19,1)
srs.SetWellKnownGeogCS('NAD83')    
ds.SetProjection( srs.ExportToWkt() )

# write the data to a single band
ds.GetRasterBand(1).WriteArray(zi)
# close
ds = None
share|improve this answer
add comment

I have also had problems with XYZ files not being recognized in gdal. I don't necessarily think that you did anything wrong. GDAL will only open gridded XYZ data files (see here). Ungridded XYZ files need to be handled in some other way. I suspect (but don't know for sure) that what SAGA has exported is an ungridded XYZ format.

In my particular use-case, I had ungridded headerless XYZ files. I worked around the problem by opening the XYZ file in a text editor and adding a header line in the format: x,y,z. I then treated the file as a csv point dataset and used gdal_rasterize to convert to a raster. I had a lot of files so I wrote a little Python script as a batch process. I hope that helps.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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