使用gdal python从geotiff读取高程

7

我正在使用GDAL加载geotiff文件。我已经成功读取了X,Y坐标,但无法读取海拔高度。

之前有人遇到过类似的情况吗?

谢谢!

1个回答

26

如果您想将所有高程值读入numpy数组中,通常会执行以下操作:

from osgeo import gdal
gdal.UseExceptions()

ds = gdal.Open('test_data.tif')
band = ds.GetRasterBand(1)
elevation = band.ReadAsArray()

print elevation.shape
print elevation

elevation将是一个2D的numpy数组。如果您想快速绘制这些值,可以使用matplotlib

import matplotlib.pyplot as plt
plt.imshow(elevation, cmap='gist_earth')
plt.show()

enter image description here

如果您想看到带有正确的x,y坐标的图形,则可以执行类似于以下内容的操作:

nrows, ncols = elevation.shape

# I'm making the assumption that the image isn't rotated/skewed/etc. 
# This is not the correct method in general, but let's ignore that for now
# If dxdy or dydx aren't 0, then this will be incorrect
x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()

x1 = x0 + dx * ncols
y1 = y0 + dy * nrows

plt.imshow(elevation, cmap='gist_earth', extent=[x0, x1, y1, y0])
plt.show()

enter image description here


2
您也可以使用 elevation = ds.ReadAsArray() 一次性读取所有波段。 - Rutger Kassies
2
Joe,请更正你的代码,将xmin改为x0以避免混淆。我不能更改它,因为编辑需要至少更改6个字符,而这不是这种情况;) - iblasi
1
@iblasi - 感谢你的发现! - Joe Kington

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接