I would like to get the coordinates X Y Z from a geotiff file.

I am using the gdal translate for the conversion.

Which output format is the more adequate for containing a 3 numpay array X Y Z ?

  • Can you clarify what you are trying to achieve? It sounds like you have a DTM. In which case you can just open the raster as an array in Numpy. this will give you an array of Z values. The x,y are then derived by the column/row in relation to the raster origin and resolution. No conversion is necessary (i.e. unless you are doing something unusual, you only need a 2-dimensional array). – MappaGnosis Jul 31 '14 at 15:35
  • I have used a gdal script before , I have managed to read X Y Z , for a a large geotiff it takes a lot of a time. I am trying to do it using the gdal translate tool – user3548574 Aug 1 '14 at 7:39

To convert GeoTIFF to XYZ you can simply use:

gdal_translate -of XYZ geotiff.tif geotiffas.xyz

From python for all tif files in the directory:

def main(my_path):
  import os
  import subprocess

  # What/Where to load and save
  source_directory = my_path
  target_directory = my_path
  # Create dir for processed files
  if not os.path.exists(target_directory):
  # Find all tif
  file_names = [files for files in os.listdir(source_directory) if files.endswith(".tif")]

  for file_name in file_names:

      # Convert tif to XYZ
      gdal_command = "gdal_translate -of XYZ \"" + source_directory + "\\" + file_name + "\" \"" + target_directory + "\\" + file_name[:-3] + "xyz" + "\""

  print("All done!")

# Call itself automatically only if run as stand alone  
if __name__ == "__main__":
  import os
  my_path = os.getcwd()

HERE is what I have managed to do

band = ds.GetRasterBand(1)

width = ds.RasterXSize
height = ds.RasterYSize

gt = ds.GetGeoTransform()
minX = gt[0]
minY = gt[3] + width*gt[4] + height*gt[5] 
maxX = gt[0] + width*gt[1] + height*gt[2]
maxY = gt[3] 

print "the domain :" , "[" ,minX,";",maxX,"]","[", minY,";",maxY ,"]"

#showing a 2D image of the topo
# plt.imshow(elevation, cmap='gist_earth',extent=[minX, maxX, minY, maxY])
# plt.show()

# elevation 2D numpy array
elevation = band.ReadAsArray()

a = gt[1]
b = gt[2]
d = gt[4]
e = gt[5]

def pixelcoord(x, y):
  """Returns coordinates X Y from pixel"""
  xp = a * x + b * y + minX
  yp = d * x + e * y + minY
  return xp, yp

#"creating X Y Z file..."
f = open("coordinates.dat", "wb")
for i in range(height):
  for j in range(width):
     pp= pixelcoord(i,j)
     print >> f, pp[0] , pp[1] , elevation[i][j]
  • Warning: pixelcoord gives the pixel corner coordinate (usually upper left). To get the centre, you need to apply a half-pixel transformation. – Mike T Jun 30 '15 at 4:59
  • Is there any reason to write this script when you can simply use gdal_translate -of XYZ geotiff.tif geotiffas.xyz ? – Miro Aug 2 '15 at 5:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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