matplotlib 从.tiff绘制3D曲面-为什么它看起来这么奇怪?

brtdzjyr  于 2023-03-30  发布在  其他
关注(0)|答案(1)|浏览(147)

我试图用matplotlib将两个.tiffs可视化为3D曲面图。结果看起来有点奇怪-有人能帮助我吗?

# import relevant moduls
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

# Specify input files
dem_old = '../data/Nisqually_1987.tif'
dem_new = '../data/Nisqually_2017.tif'

# Creating first 3D figure
# set up a figure twice as wide as it is tall
fig = plt.figure(figsize=plt.figaspect(0.5))

# Old DEM plot
ax = fig.add_subplot(1, 2, 1, projection='3d')
dem = gdal.Open(dem_old)
dem_array = dem.ReadAsArray()

lin_x = np.linspace(0,1,dem_array.shape[0],endpoint=False)
lin_y = np.linspace(0,1,dem_array.shape[1],endpoint=False)
y,x = np.meshgrid(lin_y,lin_x)
z = dem_array
surf = ax.plot_surface(x,y,z,cmap='terrain', edgecolor='none')
fig.colorbar(surf, shrink=0.5, aspect=15)
ax.set_title('3D Nisqually-Gletscher Oberfläche (1987)')
plt.xticks([])
plt.yticks([])

# New DEM 3D plot
ax = fig.add_subplot(1, 2, 2, projection='3d')
dem_2 = gdal.Open(dem_new)
dem_array_2 = dem_2.ReadAsArray()

lin_x_2 = np.linspace(0,1,dem_array_2.shape[0],endpoint=False)
lin_y_2 = np.linspace(0,1,dem_array_2.shape[1],endpoint=False)
y,x = np.meshgrid(lin_y_2,lin_x_2)
z = dem_array_2
surf = ax.plot_surface(x,y,z,cmap='terrain', edgecolor='none')
fig.colorbar(surf, shrink=0.5, aspect=10)
ax.set_title('3D Nisqually-Gletscher Oberfläche (2017)')
plt.xticks([])
plt.yticks([])

# show plot
plt.savefig ('glaciers_3D.png')
plt.show()

数据来自:https://www.sciencebase.gov/catalog/item/6140d0a3d34e1449c5d6000a

7jmck4yq

7jmck4yq1#

在检查tiff中的值后,如图所示,我发现正常值在[1000, 3000]之间,但有些值是-3.40282e+38的常数,因此只需通过以下操作将这些值更改为np.nan

dem_array = dem.ReadAsArray()
dem_array[abs(dem_array) > 1e10] = np.nan

dem_array_2 = dem_2.ReadAsArray()
dem_array_2[abs(dem_array_2) > 1e10] = np.nan

然后你会得到这个数字:

但是在2017年的图中,我发现曲面图中有一些“洞”,所以我用pandas进行了线性插值来填充“洞”:

import pandas as pd 

z = dem_array_2
df = pd.DataFrame.from_records(z)
df_padded  = df.interpolate(limit_area='inside')
z = df_padded.values

我们得到新的图形:

相关问题