netCDF4(.nc)文件读取转为tif或csv（python）

.nc文件（network Common Data Format）文件是气象上常用的数据格式，python上读取.nc使用较多的库为netCDF4这个库。基本操作参考fangzuliang的博客

"""
author: shuaijie
intro: 
date: 08/03/2020 18:47
"""
import netCDF4 as nc
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import matplotlib
matplotlib.rcParams['backend'] = 'SVG'


def main():
    var = 'CL'
    filename = r'C:\Users\13290\Desktop\soil data\{}.nc'.format(var)
    f = nc.Dataset(filename)
    all_vars_info = f.variables.items()  # 获取所有变量信息，包括名字，数据形状，以便于储存分析
    print(all_vars_info)
    # var_info = f.variables[var]  # 获取变量信息
    # 接下来的是我自己数据的特征处理
    var_data0 = f[var][0, :]  # 获取变量的数据，0代表第一层，我的是三维数据，选择第一层，每层是个二维
    var_data1 = f[var][1, :]
    var_data2 = f[var][2, :]
    var_data3 = f[var][3, :]
    var_lon = f['lon'][:]
    var_lat = f['lat'][:]
    data_array0 = pd.DataFrame(var_data0, index=var_lat, columns=var_lon)
    data_array0[data_array0 == -999] = None
    data_array0.sort_index(inplace=True, ascending=False)
    data_array1 = pd.DataFrame(var_data1, index=var_lat, columns=var_lon)
    data_array1[data_array1 == -999] = None
    data_array1.sort_index(inplace=True, ascending=False)
    data_array2 = pd.DataFrame(var_data2, index=var_lat, columns=var_lon)
    data_array2[data_array2 == -999] = None
    data_array2.sort_index(inplace=True, ascending=False)
    data_array3 = pd.DataFrame(var_data3, index=var_lat, columns=var_lon)
    data_array3[data_array3 == -999] = None
    data_array3.sort_index(inplace=True, ascending=False)
    data_array = (data_array0*4.5+data_array1*4.6+data_array2*7.5+data_array3*12.3)/28.9
    # print(data_array.head())
    # np.savetxt('SA.csv', data_array, delimiter=',')  # 可以用numpy存为.csv
    # data_array.to_excel('SA.xlsx')  # 或者存成Excel
    plt.figure(figsize=(7.56, 4.32))
    camp = sns.light_palette("darkred")
    sns.heatmap(data_array, xticklabels=False, yticklabels=False, cbar=True, robust=True,
                cbar_kws=dict(use_gridspec=False, location="right", shrink=0.4),
                cmap=camp)  # 还行吧就这，需要先画画布，再不占用画布会颜色条
    # plt.tight_layout()  # 防止保存时大小不合适
    plt.savefig(r'C:\Users\13290\Desktop\{}30.png'.format(var), dpi=500)


if __name__ == '__main__':
    main()

增强版，选取特定范围的数据，之后输出为tif或者csv

import netCDF4 as nc
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from osgeo import gdal,osr


var = 'SOM'
filename = r'C:\Users\Administrator\Desktop\soil_data\{}.nc'.format(var)
f = nc.Dataset(filename)
all_vars_info = f.variables.items()  # 获取所有变量信息，包括名字，数据形状，以便于储存分析
print(all_vars_info)

var_data0 = f[var][0, :]  # 获取变量的数据，0代表第一层，我的是三维数据，选择第一层，每层是个二维
var_data1 = f[var][1, :]
var_data2 = f[var][2, :]
var_data3 = f[var][3, :]
var_data4 = f[var][4, :]
var_data5 = f[var][5, :]
var_data6 = f[var][6, :]
var_lon = f['lon'][:]
var_lat = f['lat'][:]
data_array0 = pd.DataFrame(var_data0, index=var_lat, columns=var_lon)
data_array0[data_array0 == -999] = None
data_array0.sort_index(inplace=True, ascending=False)
data_array1 = pd.DataFrame(var_data1, index=var_lat, columns=var_lon)
data_array1[data_array1 == -999] = None
data_array1.sort_index(inplace=True, ascending=False)
data_array2 = pd.DataFrame(var_data2, index=var_lat, columns=var_lon)
data_array2[data_array2 == -999] = None
data_array2.sort_index(inplace=True, ascending=False)
data_array3 = pd.DataFrame(var_data3, index=var_lat, columns=var_lon)
data_array3[data_array3 == -999] = None
data_array3.sort_index(inplace=True, ascending=False)
data_array4 = pd.DataFrame(var_data4, index=var_lat, columns=var_lon)
data_array4[data_array4 == -999] = None
data_array4.sort_index(inplace=True, ascending=False)
data_array5 = pd.DataFrame(var_data5, index=var_lat, columns=var_lon)
data_array5[data_array5 == -999] = None
data_array5.sort_index(inplace=True, ascending=False)
data_array6 = pd.DataFrame(var_data6, index=var_lat, columns=var_lon)
data_array6[data_array6 == -999] = None
data_array6.sort_index(inplace=True, ascending=False)
# 对所有层求平均
data_array = (data_array0*4.5+data_array1*4.6+data_array2*7.5+data_array3*12.3+data_array4*20.4+
             data_array5*33.6+data_array6*55.4)/138.3
# print(data_array.head())
# 获取河南省位置的数据
left_lon = np.where(110<var_lon)[0][0]  # 左边坐标
right_lon = np.where(117>var_lon)[0][-1]  # 右边坐标
up_lat = np.where(31<var_lat[::-1])[0][-1]
low_lat = np.where(36.5>var_lat[::-1])[0][0]
# 输出为csv
data_array_henan = data_array.iloc[low_lat:up_lat, left_lon:right_lon]
data_array_henan.to_csv('henan%s.csv'%var)  # 可以用pandas存为.csv

之后可以利用输出的csv文件转为tif,也可以直接输出tif,此处以读文件为例

# 读取数据
read = pd.read_csv('henanSOM.csv')
read=read.set_index('Unnamed: 0')
# 导出为tif
var_lon = read.columns.map(float)  # 列索引为字符串，转为数字
var_lat = read.index
data_arr = np.asarray(read)
LonMin, LatMax, LonMax, LatMin = [var_lon.min(), var_lat.max(), var_lon.max(), var_lat.min()]
# 分辨率计算
N_Lat = len(var_lat)
N_Lon = len(var_lon)
Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
 
# 创建.tif文件
driver = gdal.GetDriverByName('GTiff')
out_tif_name = r'henan{}.tif'.format(var)
out_tif = driver.Create(out_tif_name, N_Lon, N_Lat, 1, gdal.GDT_Float32)  # 创建框架
 
# 设置影像的显示范围
# Lat_Res一定要是-的
geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
out_tif.SetGeoTransform(geotransform)
 
# 获取地理坐标系统信息，用于选取需要的地理坐标系统
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)  # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
out_tif.SetProjection(srs.ExportToWkt())  # 给新建图层赋予投影信息
 
# 数据写出
out_tif.GetRasterBand(1).WriteArray(data_arr)  # 将数据写入内存，此时没有写入硬盘
out_tif.FlushCache()  # 将数据写入硬盘
out_tif = None  # 注意必须关闭tif文件