python

paopao123

我用python将nc降雨数据转tif数据的时候，他一直报错数据维度有问题，我降雨变量的维度是（1，114，96），真心恳请大家帮我看看代码哪里出了问题，恳请气象学大佬给点建议

这个是他的代码呢
import os
from osgeo import gdal
import netCDF4 as nc
import numpy as np
from glob import glob
from osgeo import osr

# 导入所需库
# os 用于基本操作系统操作
# gdal 用于处理地理空间数据
# netCDF4 用于读取 netCDF 文件
# numpy 用于数值操作
# glob 用于文件路径匹配
# osr 用于处理空间参考信息

# 定义工作路径和输出路径
WorkPath = r'D:\WRFOUT2000-2020\2001\shiyan1'
OutPath = WorkPath

# 指定降雨类型
#rain = '__xarray_dataarray_variable__'

# 定义空间分辨率
SP = 0.1  # 度

# 如果输出路径不存在，则创建
if not os.path.exists(OutPath):
    os.makedirs(OutPath)

# 使用 glob 查找所有 .nc 文件
path = glob(os.path.join(WorkPath, '*.nc'))

# 对每个文件进行处理
for file in path:
    f = nc.Dataset(file)  # 打开 netCDF 文件
    #降雨的输出变量有四个
    rain1=f['RAINC'][:]
    rain2=f['RAINNC'][:]
    rain3=f['I_RAINNC'][:]*100
    rain4=f['I_RAINNC'][:]*100
    rain=rain1+rain2+rain3+rain4#wrf输出的总降水
    data = np.asarray(rain)  # 提降雨数据数据
    #data[data == 65535] = np.nan  # 将值为 65535 的数据替换为 NaN

    lon = np.array(f['XLONG'][:])  # 提取经度数据
    lat = np.array(f['XLAT'][:])
  #  lon = lon.ravel()
    #lat = lat.ravel()
    # 提取纬度数据
    LonMin, LatMax, LonMax, LatMin = lon.min(), lat.max(), lon.max(), lat.min()  # 计算范围和极值
    N_Lat = len(lat)
    N_Lon = len(lon)
    Lon_Res = SP
    Lat_Res = SP

    # 提取文件名并构建输出路径
    fname = os.path.basename(file).split('.nc')[0]
    outfile = OutPath + '/{}.tif'.format(fname)

    driver = gdal.GetDriverByName('GTiff')  # 获取 GeoTIFF 驱动
    # 创建输出 GeoTIFF 文件
    outRaster = driver.Create(outfile, N_Lon, N_Lat, 1, gdal.GDT_Float32)
    # 设置地理变换参数
    outRaster.SetGeoTransform([LonMin - Lon_Res / 2, Lon_Res, 0, LatMax + Lat_Res / 2, 0, -Lat_Res])
    sr = osr.SpatialReference()
    sr.SetWellKnownGeogCS('WGS84')  # 使用 WGS84 坐标系
    outRaster.SetProjection(sr.ExportToWkt())  # 设置投影信息
    outRaster.GetRasterBand(1).WriteArray(data)  # 将数据写入栅格波段
    print(fname + '.tif', '处理完成')

    # 释放内存
    del outRaster
    f.close(）

墨家大宝

可能你的data是三维（1，114，96）的，把你的`outRaster.GetRasterBand(1).WriteArray(data)`改成`outRaster.GetRasterBand(1).WriteArray(data[0])`试试

paopao123

墨家大宝 发表于 2024-4-18 08:48
可能你的data是三维（1，114，96）的，把你的`outRaster.GetRasterBand(1).WriteArray(data)`改成`outRaste ...

谢谢老哥，我试试

paopao123

墨家大宝 发表于 2024-4-18 08:48
可能你的data是三维（1，114，96）的，把你的`outRaster.GetRasterBand(1).WriteArray(data)`改成`outRaste ...

大佬，我找到原因了，因为我的经纬度坐标全是二维的， N_Lat = len(lat)， N_Lon = len(lon)，这两行代码构建了（1*1）的数组，而我实际数据需要一个（96*114),所以说一直报错维度不同