#導(dǎo)入庫
import netCDF4
from netCDF4 import Dataset
 
#讀取數(shù)據(jù)文件
nc_file=Dataset("/media/hsy/HSYWS/001DATA/VPD_DATA/vapor_pressure_deficit_1979.nc")
 
#輸出數(shù)據(jù)文件的兩種方式
#nc_file
print(nc_file)

輸出結(jié)果展示：

<class 'netCDF4._netCDF4.Dataset'> root group (NETCDF4 data model, file format HDF5): Conventions: CF-1.4

created_by: R, packages ncdf4 and raster (version 3.3-13)

date: 2021-10-08 13:21:50

dimensions(sizes): Longitude(1440), Latitude(721), Time(365)

variables(dimensions): int32 crs(), float64 Longitude(Longitude), float64 Latitude(Latitude), int32 Time(Time), float32 VPD(Time, Latitude, Longitude)

groups:

#所有變量讀取
print(nc_file.variables.keys())
 
#輸出結(jié)果：dict_keys(['crs', 'Longitude', 'Latitude', 'Time', 'VPD'])
 
#單個變量讀取
nc_file['Longitude']
#print(nc_file.variables['Longitude'])
 
"""<class 'netCDF4._netCDF4.Variable'>
float64 Longitude(Longitude)
    units: degrees_east
    long_name: Longitude
unlimited dimensions: 
current shape = (1440,)
filling on, default _FillValue of 9.969209968386869e+36 used"""
 
print(nc_file.variables['crs'])

結(jié)果輸出解讀：

<class 'netCDF4._netCDF4.Variable'> #文件數(shù)據(jù)類型

int32 crs()

proj4: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs #proj4坐標系參數(shù)，詳情請見：Quick start — PROJ 9.1.0 documentation

unlimited dimensions:

current shape = () filling on, default _FillValue of -2147483647 used

#單個變量的所有屬性名稱
print(nc_file.variables['VPD'].ncattrs())
#['_FillValue', 'long_name', 'grid_mapping', 'proj4', 'min', 'max']
print(nc_file.variables['VPD'].proj4)#proj4坐標系
print(nc_file.variables['VPD'].grid_mapping)#給定坐標變量與真實經(jīng)緯度坐標之間的映射關(guān)系：crs
print(nc_file.variables['VPD']._FillValue)#填充值或空值
#讀取變量的維度
print(nc_file['VPD'].shape)
#(365, 721, 1440) （#time, latitude, longitude) 格點分辨率為：天*0.25*0.25度。
#讀取變量值
VPD=nc_file.variables['VPD'][:]
print(VPD)#讀取結(jié)果含有全部數(shù)值。
#print(nc_file['VPD'])#輸出結(jié)果不完整

繪圖

1、利用matplotlib繪圖

import matplotlib.pyplot as plt
plt.contourf(long, lat, VPD[10, :, :] )
plt.colorbar(label="VPD", orientation="horizontal")
plt.show()

2、利用basemap繪圖

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
lon0 = long.mean()
lat0 = lat.mean()
# 設(shè)置投影方式：cyl為圓柱投影、還可設(shè)置merc為mercator投影 llcrnrlat為起始lat；urcrnrlat為終止lat
# m = Basemap(projection='merc', llcrnrlat=lat[0], urcrnrlat=lat[-1], \
#              llcrnrlon=lon[0], urcrnrlon=lon[-1], ax=ax1)
# 參數(shù) "resolution" 用于控制地圖面積邊緣的精細程度，有'l'和'h'兩種取值
m = Basemap(lat_0=lat0, lon_0=lon0,projection='cyl',resolution='l')
# 繪制等經(jīng)緯度線 緯度每隔20度畫一條線，且標注經(jīng)緯度
m.drawparallels(np.arange(-90., 91., 20.), labels=[1, 0, 0, 0], fontsize=10)
m.drawmeridians(np.arange(-180., 181., 40.), labels=[0, 0, 0, 1], fontsize=10)
m.drawcoastlines()# 繪制海岸線
lon, lat = np.meshgrid(long, lat)
xi, yi = m(lon, lat)
# cmap是圖形顏色，還可選‘jet'、‘spring'、‘winter'、'summer'、'autumn'
cs = m.contourf(xi, yi, VPD[10],  cmap='summer')
# pad指位置，
cbar = m.colorbar(cs, location='bottom', pad="10%",format='%.1f')
font1 = {'family': 'DejaVu Sans', 'weight': 'normal', 'size': 16}
plt.title('VPD', font1)
plt.show()

3、利用cartopy繪圖

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(15, 7))  
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': proj})  
# 或者 ax = fig.add_subplot(111,proj=proj)
lon1 = nc_file.variables['Longitude'][:]
lat1 = nc_file.variables['Latitude'][:]
print(lon1.shape, lat1.shape)
ax.contourf(lon, lat, VPD[100])
ax.coastlines(resolution = '10m')
#添加格網(wǎng)
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
# 設(shè)置 gridlines 和 ticklabels 
gl = ax.gridlines(draw_labels = True, linewidth = 1.5)
gl.xlabels_top = False 
gl.xlines = True
gl.xformatter = LONGITUDE_FORMATTER
gl.ylabels_right = False 
gl.ylines = True
gl.yformatter = LATITUDE_FORMATTER
plt.show()