python批量解析CLDAS-V2.0的.nc数据文件

原创已于 2023-04-18 18:23:00 修改 · 5.1k 阅读

50 ·

CC 4.0 BY-SA版权

文章标签：

#NetCDF #气象数据 #python

于 2020-12-16 15:23:00 首次发布

地理数据专栏收录该内容

10 篇文章

订阅专栏

这篇博客介绍了如何使用Python的NetCDF4和gdal库将CLDAS-V2.0的.nc数据文件（以四川省为例）转换为栅格图像。首先通过ArcMap的多维工具预览数据，然后通过编写Python脚本读取.nc文件，利用NetCDF4获取变量，并用gdal创建栅格数据，设置地理信息和投影，最后将数据写入.tif文件。脚本适用于批量处理多个.nc文件，并可在ArcMap中正常显示结果。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

处理CLDAS-V2.0的.nc数据文件——四川省为例

CLDAS-V2.0是中国气象局陆面数据同化系统实时产品数据集，可以从中国气象数据网下载在这里插入图片描述
以四川省为例，下载了一段时间的数据，数据格式为.nc格式（NETCDF DATA），如下

可以使用ArcMap的多维工具的“创建NetCDF栅格图层”的工具来打开，叠加一个全国的省界，可以看到这个数据是四川的某个数据
在这里插入图片描述
现在使用NetCDF4和gdal来处理这个文件夹中的459个.nc文件，转为栅格图像

使用NetCDF4读取一个nc文件，可以看到是这样的

其中variables是它的值，有经纬度（坐标）和一个TAIR的值（灰度）
下面pip好NetCDF4和gdal两个包，开始写一个处理的py文件

import netCDF4
from osgeo import gdal
import numpy as np
import os

# 读取一个nc文件转换为栅格
def read_nc_to_tif(nc_path, out_path):
    # 读取nc数据保存为字典
    dataset = netCDF4.Dataset(nc_path)  # 打开一个nc文件
    keys = dataset.variables.keys()  # 获取它的值的列表
    dicts = {}
    for key in keys:
        dicts[key] = np.array(dataset[key])  # 保存为字典
    # 创建栅格数据
    driver = gdal.GetDriverByName("GTiff")  # 新建一个空的tif
    _, file_name = os.path.split(nc_path)  # 得到nc文件的文件名
    outdata = driver.Create(
        os.path.join(out_path, (file_name.split('.')[0]+'.tif')),  # 要保存的tif文件名
        dicts['LON'].shape[0],  # x的size
        dicts['LAT'].shape[0],  # y的size
        len(list(dicts.keys()))-2,  # 波段数（除去经纬度其他的都是波段）
        gdal.GDT_Float32  # 类型，这里的这些文件都是float32的，其他的根据实际情况改
    )
    # 计算转换矩阵
    pos_x = dicts['LON'].min()
    pos_y = dicts['LAT'].min()
    pix_x = (dicts['LON'].max() - pos_x) / dicts['LON'].shape[0]
    pix_y = (dicts['LAT'].max() - pos_y) / dicts['LAT'].shape[0]
    outdata.SetGeoTransform((pos_x, pix_x, 0, pos_y, 0, pix_y))
    # 设置投影为WGS84
    proj_str = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG",4326]]'
    outdata.SetProjection(proj_str)
    # 投影和转换根据实际情况
    idx = 0  # 波段计数
    for ikey in dicts.keys():
        if ikey != 'LON' and ikey != 'LAT':  # 经纬度不用于产生新的波段
            outband = outdata.GetRasterBand(idx+1)  # 波段加1
            # print(str(idx), ikey)
            outband.WriteArray(dicts[ikey].astype('float32'))  # 写入数据
            idx += 1
    outdata.FlushCache()
    outdata = None
    print(nc_path, ' Finished')
    
# 主程序
if __name__ == '__main__':
    # 输入nc文件文件夹
    data_path = 'E:/XXXXXXX/NC_20201216_0950'
    # 输出栅格文件夹
    output_path = 'E:/XXXXXXX/results'
    # 处理流程
    ncs = os.listdir(data_path)
    for nc in ncs:
        nc_path = os.path.join(data_path, nc)
        read_nc_to_tif(nc_path, output_path)
    print('All Done!')