请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 283|回复: 3

python

[复制链接]

新浪微博达人勋

发表于 2024-4-17 22:59:58 | 显示全部楼层 |阅读模式
5金钱
我用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()

]397}(G]MO%NXTHNQ9{B9W0.png

nc2geotif.py

2.3 KB, 下载次数: 0

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2024-4-18 08:48:21 | 显示全部楼层
可能你的data是三维(1,114,96)的,把你的`outRaster.GetRasterBand(1).WriteArray(data)`改成`outRaster.GetRasterBand(1).WriteArray(data[0])`试试
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-4-18 09:13:44 | 显示全部楼层
墨家大宝 发表于 2024-4-18 08:48
可能你的data是三维(1,114,96)的,把你的`outRaster.GetRasterBand(1).WriteArray(data)`改成`outRaste ...

谢谢老哥,我试试
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2024-4-19 11:42:55 | 显示全部楼层
墨家大宝 发表于 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),所以说一直报错维度不同
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表