处理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!')
接下来运行脚本,即可在results文件夹下得到tif文件,但是看起来是一片空白
没事,因为是float32的类型,拖入ArcMap中,没有问题,和直接导入的完全重合,数值也一样,黑白还是彩色去显示设置设置一下即可