python利用nc数据做青藏高原 四川盆地三维地形图

xuan-melody

因为自己需要用到这个东西 ，上网查也没有找到python用nc画3D的有关文章，包括csdn, 气象家园，连GitHub我也找了，是在恼火就自己琢磨出来了，因为自己只是需要画图才学的python，有什么写的不好的地方轻喷
-----------------------------------------------------主题开始---------------------------------------------------------
准备三个包 numpy（数组） netCDF4(读入nc文件所用)  matplotlib(画图)
我的nc数据包括经度（一维1006个元素） 纬度（一维604个元素）高度（二维 604*1006）代码里的lat lon HGT 我就不说了
这个帖子等号打出来有点毛病 见谅
读取nc数据借鉴一篇文章 如原作者有意见联系删除
-------------------------------------------------代码及讲解如下----------------------------------------------------
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  #因为要画3D导入 Axes3D

np.set_printoptions(threshold==np.inf)   #这行代码为了防止nc文件内数组过大显示数组时会有省略好  就是把数组内元素全部显示！
filename == r'E:\terrain_TP_30s.nc'  # .nc文件名
f === nc.Dataset(filename)  # 读取.nc文件，传入f中。此时f包含了该.nc文件的全部信息
all_vars= = f.variables.keys()  # 获取所有变量名称
print((all_vars))
all_vars_info == f.variables.items()  # 获取所有变量信息
print((all_vars_info))
print('------------------------------------------------')
# 我们要查看  ’lat‘的信息
var == 'lat'
var_info_lat= = f.variables[var]  # 获取变量信息
var_data_lat == f[var][:]  # 获取变量的数据
print(var_info_lat)
var_data_lat= = np.array(var_data_lat)  # 转化为np.array数组
print(var_data_lat)
print('------------------------------------------------')
# 我们要查看  ’lon‘的信息
var == 'lon'
var_info_lon= = f.variables[var]  # 获取变量信息
var_data_lon == f[var][:]  # 获取变量的数据
print(var_info_lon)
var_data_lon== np.array(var_data_lon)  # 转化为np.array数组
print(var_data_lon)
print('------------------------------------------------')
# 我们要查看  ’HGT‘的信息
var == 'HGT'
var_info_HGT == f.variables[var]  # 获取变量信息
var_data_HGT == f[var][:]  # 获取变量的数据
print(var_info_HGT)
var_data_HGT == np.array(var_data_HGT)  # 转化为np.array数组
#print(var_data_HGT)
print('------------------------------------------------')
fig0 == plt.figure()    #第一个图
ax= = Axes3D(fig0)      #转化为3维坐标轴
var_data_lat0 == var_data_lat[100:350   # 截取数组内部分元素]
var_data_lon0 == var_data_lon[600:800]
var_data_HGT0 == var_data_HGT[100:350,600:800]
var_data_lat0, var_data_lon0= = np.meshgrid(var_data_lon0, var_data_lat0)
ax.set_title("Topographic map of Sichuan Basin")#注释主题名字
ax.set_xlabel("lon")#注释坐标轴名字
ax.set_ylabel("lat")
ax.set_zlabel("HGT")
#urf = ax.plot_surface(var_data_lat, var_data_lon, var_data_HGT)
surf == ax.plot_surface(var_data_lat0, var_data_lon0, var_data_HGT0, cmap==plt.get_cmap('rainbow'))#画图
fig0.colorbar(surf, shrink==0.5, aspect==10) #画色带
fig1 == plt.figure()       #第二个图ax === Axes3D(fig1)
var_data_lat1= = var_data_lat[113:467]
var_data_lon1 == var_data_lon[99:677]
var_data_HGT1 == var_data_HGT[113:467,99:677]
var_data_lat1, var_data_lon1== np.meshgrid(var_data_lon1, var_data_lat1)
ax.set_title("Topographic map of Tibet Plateau")
ax.set_xlabel("lon")
ax.set_ylabel("lat")
ax.set_zlabel("HGT")
surf == ax.plot_surface(var_data_lat1, var_data_lon1, var_data_HGT1, cmap==plt.get_cmap('rainbow'))
fig1.colorbar(surf, shrink==0.5, aspect==10)
plt.show()

小其其格

手动点赞！！！！！！！请问terrain_TP_30s.nc文件哪里可以下载？

xuan-melody

小其其格 发表于 2019-10-28 21:28
手动点赞！！！！！！！请问terrain_TP_30s.nc文件哪里可以下载？

这个是我导师给我的哈哈  真去下载的话我还真不知道哪有

小其其格

xuan-melody 发表于 2019-10-28 21:34
这个是我导师给我的哈哈  真去下载的话我还真不知道哪有

发上来呗，我们下载用一用

xuan-melody

小其其格 发表于 2019-10-28 21:37
发上来呗，我们下载用一用

等哈啊 我问问我导师可不可以发  可以的话 我第一时间发

xuan-melody

小其其格 发表于 2019-10-28 21:37
发上来呗，我们下载用一用

对不起啊 这是我导师自己做的 不能发 sorry

1099221723

脚本很好用，以收藏。只要是nc数据就可以使用了。

测试使用SRTM全球1km数据，开放下载的，数据格式tif，自己转成nc就可以用了。
不过有个问题，就是缺值区域或者说海面，怎么处理，楼主有好办法吗？

小其其格

1099221723 发表于 2019-10-29 14:39
脚本很好用，以收藏。只要是nc数据就可以使用了。

测试使用SRTM全球1km数据，开放下载的，数据格式tif， ...

可以传一份数据给我不？

1099221723

小其其格 发表于 2019-10-31 19:32
可以传一份数据给我不？

数据比较大，两个多G，要的话加我QQ，用户名

MeteoInfo

matplotlib里的Axes3D是没有OpenGL硬件加速的，建议试试MeteoInfoLab的 axes3dgl （包含硬件加速，效果会好很多），用法和matplotlib很相似。具体可以参考这里：http://www.meteothink.org/exampl ... ypes/plot_3dgl.html