python gdal nc数据转tif 包括如何获取变量及变量属性

news2024/11/23 14:38:21

文章目录

1 gdal nc转换为tif
2 利用netCDF4获取变量及变量属性
3 进行gregorian时间类型的转换
4 总结

1 gdal nc转换为tif

地理变换是从图像坐标空间（行、列）（也称为（像素、线）到地理参考坐标空间（投影或地理坐标）的仿射变换。

总共包含6个元素

ExportToWkt 是 GDAL（Geospatial Data Abstraction Library）库中的一个函数，用于将空间参考信息导出为 Well-Known Text（WKT）格式。

Well-Known Text（WKT）是一种文本表示形式，用于描述空间参考系统（Spatial Reference System）的信息，包括坐标系、投影和地理转换参数等。WKT 是一种通用的标准，可以在不同的 GIS（地理信息系统）和地理数据处理软件之间共享和交换空间参考信息。
这里以ERA5数据为例（按月获取小时温度数据）

def nc_to_tif(data, lons, lats, timestamp):
    # 影像的左上角和右下角坐标
    LonMin, LatMax, LonMax, LatMin = [lons.min(), lats.max(), lons.max(), lats.min()]
    lon_number = len(lons)
    lat_number = len(lats)
    # 1 创建tif文件
    driver = gdal.GetDriverByName('GTiff')
    out_tif_name = output_directory + timestamp + '_t2m.tif'
    out_tif = driver.Create(out_tif_name, lon_number-1, lat_number-1, 1, gdal.GDT_Float32)

    # 2 设置影像显示范围[建立仿射变换]
    # 分辨率与原始数据的空间分辨率需要保持一致
    Lon_Res = 0.1
    Lat_Res = 0.1
    geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
    out_tif.SetGeoTransform(geotransform)
    # 3 将地理经纬度坐标转换为投影坐标
    projection = gdal.osr.SpatialReference()
    # 4326指'WGS84'
    projection.ImportFromEPSG(4326)
    # 4 为tif数据赋予投影信息
    # 将空间参考以WKT格式输出
    out_tif.SetProjection(projection.ExportToWkt())
    # 5 将数据写入栅格
    out_tif.GetRasterBand(1).WriteArray(data)
    out_tif.GetRasterBand(1).SetNoDataValue(-9999)
    del out_tif

2 利用netCDF4获取变量及变量属性

from netCDF4 import Dataset

# 打开nc文件
nc_file = Dataset("test.nc", "w", format="NETCDF4")
# 获取NC文件的维度信息
print("\n维度信息:")
    for dimname, dim in nc_file.dimensions.items():
    	print(f"{dimname}: {len(dim)}")

# 获取变量信息以及变量的属性信息
print("\n变量信息:")
    for varname, var in nc_file.variables.items():
        print(f"{varname}:")
        print(f"数据类型: {var.dtype}")
        print(f"维度: {var.dimensions}")
        print(f"属性:")
        for attr in var.ncattrs():
            print(f"    {attr}: {getattr(var, attr)}")


# 根据特定的变量名获取变量
data = nc_file.variables['variable_name'][:]
# 根据特定的变量名和特定的属性名获取[对应的属性]
scale_factor = getattr(nc_file.variables['variable_name'], 'scale_factor')

3 进行gregorian时间类型的转换

公历（Gregorian calendar），也称为格里高利历、格列高利历、格雷戈里曆、西曆或新曆，是今天世界上使用最广泛的日历系统。它是由教皇格里高利十三世在1582年颁布的，以修正尤利乌斯日历的缺陷。公历采用365天5小时48分46秒的长度作为一年，并在必要时通过闰年来调整。与其他日历系统相比，公历更为准确，被广泛用于全球大部分国家和文化中。

简而言之，当提到格里高利历（Gregorian calendar）或公历时，通常是指今天全球通用的日期系统，用于标记时间、日期和年份。

进行nc时间维度的改变

利用nc.num2date()将[时间的数值数据]转换为[datetime]对象

 time_new = nc.num2date(nc_file.variables['time'][:], 'hours since 1900-01-01 00:00:00.0').data

4 总结

对nc数据处理时，需要先判断是否是_FillValue
很重要: 使用netCDF4进行数据处理时，会自动进行【data×scale_factor+add_offset】的操作，自己不需要再进行处理
这里遇到了一个问题：按照左上角点的经纬度坐标、行数和列数，最终得到的tif数据的右下角点【经度总是多一份，纬度总是小一份】
为了避免错误，我用IDL也测试了一下，也是同样的经纬度范围，这个问题先放在这