1

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 ?

2
  • 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). Commented Jul 31, 2014 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 Commented Aug 1, 2014 at 7:39

2 Answers 2

1

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):
    os.makedirs(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" + "\""
      subprocess.call(gdal_command)

  print("All done!")

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

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]
2
  • Warning: pixelcoord gives the pixel corner coordinate (usually upper left). To get the centre, you need to apply a half-pixel transformation.
    – Mike T
    Commented Jun 30, 2015 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
    Commented Aug 2, 2015 at 5:55

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.