python将nc文件(NetCDF)转成tif

原创

echohye 2023-11-24 09:37:41 博主文章分类：Python学习 ©著作权

文章标签 数据 sed 倒序 文章分类 HarmonyOS 后端开发

©著作权归作者所有：来自51CTO博客作者echohye的原创作品，请联系作者获取转载授权，否则将追究法律责任

NetCDF介绍

NetCDF（网络公用数据格式）是一种用来存储温度、湿度、气压、风速和风向等多维科学数据（变量）的文件格式。
https://zhuanlan.zhihu.com/p/600050278

https://www.osgeo.cn/gdal/drivers/vector/netcdf.html

nc转tif

安装netcdf 包

pip install netCDF4

读取netcdf数据

import netCDF4 as nc

data = r'data/asi-AMSR2-s6250-20231014-v5.4.nc'
NC_DS = nc.Dataset(data)
print(NC_DS.variables.keys())
print(NC_DS.variables) # 了解变量的基本信息

点击查看打印输出

dict_keys(['polar_stereographic', 'x', 'y', 'z'])

{'polar_stereographic': <class 'netCDF4._netCDF4.Variable'>
|S1 polar_stereographic()
    grid_mapping_name: polar_stereographic
    straight_vertical_longitude_from_pole: 0.0
    false_easting: 0.0
    false_northing: 0.0
    latitude_of_projection_origin: -90.0
    standard_parallel: -70.0
    long_name: CRS definition
    longitude_of_prime_meridian: 0.0
    semi_major_axis: 6378273.0
    inverse_flattening: 298.279411123064
    spatial_ref: PROJCS["NSIDC Sea Ice Polar Stereographic South",GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid",DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",SPHEROID["Hughes 1980",6378273,298.279411123064,AUTHORITY["EPSG","7058"]],AUTHORITY["EPSG","6054"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4054"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",-70],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3412"]]
    GeoTransform: -3950000 6250 0 4350000 0 -6250 
unlimited dimensions: 
current shape = ()
filling on, default _FillValue of   used, 'x': <class 'netCDF4._netCDF4.Variable'>
float64 x(x)
    standard_name: projection_x_coordinate
    long_name: x coordinate of projection
    units: m
unlimited dimensions: 
current shape = (1264,)
filling on, default _FillValue of 9.969209968386869e+36 used, 'y': <class 'netCDF4._netCDF4.Variable'>
float64 y(y)
    standard_name: projection_y_coordinate
    long_name: y coordinate of projection
    units: m
unlimited dimensions: 
current shape = (1328,)
filling on, default _FillValue of 9.969209968386869e+36 used, 'z': <class 'netCDF4._netCDF4.Variable'>
float32 z(y, x)
    long_name: z
    _FillValue: nan
    actual_range: [  0 100]
    grid_mapping: polar_stereographic
unlimited dimensions: 
current shape = (1328, 1264)
filling on}

从输出上我们可以知道数据具有4个属性'polar_stereographic', 'x', 'y', 'z'。

处理netcdf为tiff

python将nc文件(NetCDF)转成tif_数据

# -*- coding: UTF-8 -*-
import os
import netCDF4 as nc
import numpy as np
from osgeo import gdal, osr
import glob

os.environ['PROJ_LIB'] = r'D:\Software\Development\release-1930-x64-gdal-3-6-2-mapserver-8-0-0\bin\proj9\share'
os.environ['GDAL_DATA'] = r'D:\APP\Develop\PostgreSQL\14\gdal-data'


def nc2tif(data, Output_folder, i):
    pre_data = nc.Dataset(data)  # 利用.Dataset()读取nc数据
    print(pre_data.variables.keys())
    Lat_data = pre_data.variables['y'][:]
    Lon_data = pre_data.variables['x'][:]
    pre_arr = np.asarray(pre_data.variables['z'])  # 属性变量名
    print(Lat_data, Lon_data)

    # 影像的左上角&右下角坐标
    Lonmin, Latmax, Lonmax, Latmin = [Lon_data.min(), Lat_data.max(), Lon_data.max(), Lat_data.min()]

    # 分辨率计算
    Num_lat = len(Lat_data)
    Num_lon = len(Lon_data)
    Lat_res = (Latmax - Latmin) / (float(Num_lat) - 1)
    Lon_res = (Lonmax - Lonmin) / (float(Num_lon) - 1)
    # 创建tif文件
    driver = gdal.GetDriverByName('GTiff')
    out_tif_name = Output_folder + '\\' + data.split('\\')[-1].split('.')[0] + '_' + str(i + 1) + '.tif'
    out_tif = driver.Create(out_tif_name, Num_lon, Num_lat, 1, gdal.GDT_Float32)
    # 设置影像的显示范围
    # Lat_re前需要添加负号
    geotransform = (Lonmin, Lon_res, 0.0, Latmax, 0.0, -Lat_res)
    out_tif.SetGeoTransform(geotransform)
    # 定义投影
    prj = osr.SpatialReference()
    prj.ImportFromEPSG(3976)
    out_tif.SetProjection(prj.ExportToWkt())

    if Lat_data[0] <= Lat_data[-1]:  # 如果维度上的第一个值小于等于最后的的一个值就认为是倒序，就得进行数组的倒序排列，否则就是正向，不用倒序排列
        print(f'{data}行数据是倒的，现在进行矫正............')
        pre_arr = pre_arr[::-1]
        print('矫正完成...........')
    else:
        pass
    # 数据导出
    out_tif.GetRasterBand(1).WriteArray(pre_arr)  # 将数据写入内存
    out_tif.FlushCache()  # 将数据写入到硬盘
    out_tif = None  # 关闭tif文件


def main():
    Input_folder = r'data'
    Output_folder = r'data'
    # 读取所有数据
    data_list = glob.glob(os.path.join(Input_folder + '\*.nc'))
    print(os.getcwd())
    print(data_list)
    for i in range(len(data_list)):
        data = data_list[i]
        nc2tif(data, Output_folder, i)
        print(data + '转tif成功')


if __name__ == '__main__':
    main()

python将nc文件转成tif 详解

python处理NetCDF（.nc）为tiff

利用python将nc批量转为tif（可执行的通用代码，小白都能运行！）