【发布时间】:2014-09-17 08:42:18
【问题描述】:
我正在使用 GDAL 加载一个 geotiff 文件。我设法读取了坐标 X,Y,但没有读取海拔。
以前有没有人处理过类似的案例?
问候,
【问题讨论】:
我正在使用 GDAL 加载一个 geotiff 文件。我设法读取了坐标 X,Y,但没有读取海拔。
以前有没有人处理过类似的案例?
问候,
【问题讨论】:
如果您想将所有高程值读取到一个 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()
如果您想查看具有正确* 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()
【讨论】:
elevation = ds.ReadAsArray() 一次读取所有频段。
xmin 更改为 x0 以避免混淆。我无法更改它,因为编辑需要至少 6 个字符更改,但事实并非如此 ;)