GLDAS - nc转tif（多个变量提取）

同时提取多个变量只需将以下代码串联即可。

1、导包

import os,sys
import glob
import numpy as np
import netCDF4 as nc
from osgeo import gdal, osr
from pathlib import Path

2、建立相应的文件夹

#只需要修改输入和输出文件目录即可
Input_folder = ('H:\\GLDAS\\origin\\2014')
Output_folder = ('H:\\GLDAS\\mid_data\\2014')

os.makedirs(Output_folder + '\\Tair\\')
os.makedirs(Output_folder + '\\Psur\\')
os.makedirs(Output_folder + '\\Qair\\')
os.makedirs(Output_folder + '\\SWdown\\')
os.makedirs(Output_folder + '\\Wind\\')

Tair_output_folder = os.path.join(Output_folder,'Tair')
Psur_Output_folder = os.path.join(Output_folder,'Psur')
Qair_Output_folder = os.path.join(Output_folder,'Qair')
SWdown_Output_folder = os.path.join(Output_folder,'SWdown')
Wind_Output_folder = os.path.join(Output_folder,'Wind')

3、提取变量

该部分主要参考：https://www.jianshu.com/p/50868ea0fc25

提取更多的变量只需要修改下面函数中的个别参数即可

def nc2tiff(nc_url, var='Psurf_f_inst', save_name='./nc2tiff.tif'):
 """
    :param nc_url:  netcdf文件的路径
    :param var:   读取nc文件的变量名称
    :param save_name:  输出tiff文件的保存路径
    :return:
    """
 print(nc_url)
 nc_fp = nc.Dataset(nc_url)
 # 查看属性
 print(nc_fp)
 # # 查看变量
 print(nc_fp.variables)
 # 获取var数据集
 Psurf = nc_fp.variables[var][:] 
 Psurf = np.array(Psurf)
 Psurf = np.squeeze(Psurf)
 # missing value 处理
 fillvalue = nc_fp.variables[var].missing_value
 Psurf[Psurf == fillvalue] = np.nan
 #进行
 Psurf = Psurf / 100 #此位置决定输出数值的大小，是否直接提取还是进行计算
 array2raster(Psurf, 0.25, 0.25, [-180, -60], save_name)

# 数据矩阵转tiff文件
def array2raster(array, xsize, ysize, rasterOrigin, newTiffName):
 """
     array: 计算后的栅格数据
     xsize: x方向像元大小
     ysize: y方向像元大小
     rasterOrigin: 原始栅格数据Extent
     newRasterfn: 输出tif路径
    """
 cols = array.shape[1]  # 矩阵列数
 rows = array.shape[0]  # 矩阵行数
 originX = rasterOrigin[0]  # 起始像元经度
 originY = rasterOrigin[1]  # 起始像元纬度
 driver = gdal.GetDriverByName('GTiff')
 outRaster = driver.Create(newTiffName, cols, rows, 1, gdal.GDT_Float32)
 # 设置仿射矩阵参数（左上角X坐标; X方向分辨率; 旋转角度，如果图像北方朝上，该值为0; 左上角Y坐标; 旋转角度; y方向分辨率）
 outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
 # 获取数据集第一个波段，是从1开始，不是从0开始
 outband = outRaster.GetRasterBand(1)
 outband.WriteArray(array)
 spatialRef = osr.SpatialReference()
 # 代码4326表示WGS84坐标
 spatialRef.ImportFromEPSG(4326)
 outRaster.SetProjection(spatialRef.ExportToWkt())
 outband.FlushCache()
 outRaster.FlushCache()

# Input_folder = ('H:\\GLDAS\\origin\\2013\\222')
# Output_folder = ('H:\\GLDAS\\mid_data\\2013')
# Psur_Output_folder = os.path.join(Output_folder,'Psur')

os.chdir( Psur_Output_folder ) #设置当前的工作路径
data_list = glob.glob(Input_folder + '\\*.nc4')
for i in range(len(data_list)):
 netcdf_url = data_list[i]
 nc2tiff(netcdf_url, "Psurf_f_inst", save_name = 'GLDAS_3H_' + netcdf_url[-21:-13] + netcdf_url[-12:-10] + '_Psurf_f_inst_mb.tif')
 pass