代码收藏家技术教程 2024-11-27

【Python实例】Python读取并绘制nc数据

准备：安装netCDF库等

读取nc数据相关信息

绘制图形

利用basemap绘图

参考

准备：安装netCDF库等

以【1960-2020年中国1km分辨率月降水数据集】中2020年降水为例。

先在Panopoly中查看数据属性，如下：

右键绘制图形，如下：

读取nc数据相关信息

Python代码如下：

import netCDF4
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np

# 打开下载的.nc文件
ncfile = Dataset("D:/0 DataBase/CN_Prec_MonthlyMean_1km_2020.nc/CN_Prec_MonthlyMean_1km_2020.nc")
print(ncfile.variables.keys())

# 读取变量
lat = ncfile.variables['latitude'][:]
long = ncfile.variables['longitude'][:]
[X, Y] = np.meshgrid(lat, long)    # 合成矩阵，将向量latitude和longitude定义的区域转换成矩阵X和Y
pre = ncfile.variables['pre'][:]
time = ncfile.variables['time'][:]

# 绘制图形
plt.contourf(Y, X,pre[0])
plt.colorbar(label="MonthlyMean Precipitation", orientation="horizontal")
plt.show()

# 关闭nc文件
ncfile.close()

输出属性如下：

绘制图形

利用basemap绘图

Python代码如下：

#利用basemap绘图
import netCDF4
from mpl_toolkits import basemap
from netCDF4 import Dataset
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

# 打开下载的.nc文件
ncfile = Dataset("D:/0 DataBase/CN_Prec_MonthlyMean_1km_2020.nc/CN_Prec_MonthlyMean_1km_2020.nc")

# 读取变量
lat = ncfile.variables['latitude'][:]
long = ncfile.variables['longitude'][:]
[X, Y] = np.meshgrid(lat, long)    # 合成矩阵，将向量latitude和longitude定义的区域转换成矩阵X和Y
pre = ncfile.variables['pre'][:]
time = ncfile.variables['time'][:]

# 创建一个图形和一个轴
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)

# 设定标题
titleName = "MonthlyMean Precipitation"

# 初始化Basemap对象，设置地图参数
# 设置投影方式：cyl为圆柱投影、还可设置merc为mercator投影,lcc兰伯特投影等。
# 参数 "resolution" 用于控制地图面积边缘的精细程度，有 'crude', 'low' and'intermediate' resolution datasets are installed by default.
# 如果需要'high' or 'full' resolution dataset,需要`conda install -c conda-forge basemap-data-hires`.
m = Basemap(lat_1=30, lat_2=62, lon_0=103, projection='lcc', resolution='i',
            llcrnrlon=80, llcrnrlat=10, urcrnrlon=145, urcrnrlat=55)

# 设置经纬度
# 绘制纬度坐标，纬度每隔10度画一条线，且左侧和右侧标注纬度。
m.drawparallels(np.arange(10., 55., 10.), labels=[1, 1, 0, 0], fontsize=10)
# 绘制经度坐标，经度每隔10度画一条线，且上侧和下侧标注经度。
m.drawmeridians(np.arange(70., 160., 10.), labels=[0, 0, 1, 1], fontsize=10)

# 形状文件的路径
#shp_path = "D:/0 DataBase/0 GBA Database/Boundary_Projected"
# 添加中国地图图层
#m.readshapefile(shp_path,'Boundary_Projected', drawbounds=True)

# 添加大陆地图图层
m.fillcontinents(color="coral", lake_color="blue")
# 设定海洋的颜色
m.drawmapboundary(fill_color="aqua")

# 在指定位置绘制坐标点
lon_me = long.mean()
lat_me = lat.mean()
lon_do, lat_do = np.meshgrid(lon_me, lat_me)  # 生成网格点坐标矩阵

# 转换经纬度到Basemap坐标系
xi, yi = m(lon_do, lat_do)
# print(xi,yi)
m.plot(xi, yi, "ro")  # 绘制红色点
# 在指定位置标注文字
plt.text(xi, yi, "A1(%5.1fW, %3.1fN)" % (lon_me, lat_me))

# 关闭.nc文件
ncfile.close()

# 显示图形
plt.show()

图形绘制如下：