使用python进行nc转tif的3种情况解决

目录

• 前言
• 一、nc4_to_tif（多时间序列）
• 二、nc_to_tif（多时间序列）
• 三、nc_to_tif（单时间序列）
• 总结

前言

本人是位大二在读在校学生，专业为地理信息科学，因跟老师一起做项目，所以有幸接触nc数据转换为tif数据，因为在这件事情上也遇过不少坑，也花了不少时间，所以想在这里将自己的心得与学习的经验一起分享给大家，希望大家能获得帮助，同时我也会很开心的！！如果哪里说的有问题或是代码能再改进的话，希望大家能不吝赐教，评论区欢迎大家哦！！！！！！

一、nc4_to_tif（多时间序列）

话不多说，直接上代码(代码里面的注释也是挺详细的，所以就不赘述了)：

代码如下（示例）：

# -*- coding: utf-8 -*-
import numpy as np
import netCDF4 as nc
from osgeo import gdal,osr,ogr
import glob
import os
from zipfile import ZipFile
import shutil
band_name = ''
def NC_to_tiffs(data,out_path):

    '''
    这个函数里面有些地方还是可能需要更改
    比如：coord(坐标系)
    '''
    coord = 4326            #坐标系，["EPSG","4326"]，默认为4326
    nc_data_obj = nc.Dataset(data)
    #print(nc_data_obj,type(nc_data_obj))              # 了解nc的数据类型，<class 'netCDF4._netCDF4.Dataset'>
    #print(nc_data_obj.variables)                      #了解nc数据的基本信息
    key=list(nc_data_obj.variables.keys())            #获取时间，经度，纬度，波段的名称信息，这些可能是不一样的
    print('基础属性为-----',key)
    lon_size = [i for i,x in enumerate(key) if (x.find('lon'.upper())!=-1 or x.find('lon'.lower())!=-1)][0]   #模糊查找属于经度的名称
    lat_size = [i for i,x in enumerate(key) if (x.find('lat'.upper())!=-1 or x.find('lat'.lower())!=-1)][0]   #模糊查找属于纬度的名称
    global band_name
    if band_name == '':
        band_name = input("请输入您想要输出的波段的名字（您可以从'基础属性中得来',不用加上引号）———————:")      #这里是从用户那传入一个波段的字符串，因为nc4的数据比nc要复杂，所以要让用户确定波段的名字
    band_size = [i for i,x in enumerate(key) if (x.find(str(band_name).upper())!=-1 or x.find(str(band_name).lower())!=-1)][0]
    key_band = key[band_size]            #获取波段的名称
    key_lon = key[lon_size]              #获取经度的名称
    key_lat = key[lat_size]              #获取纬度的名称
    Lon = nc_data_obj.variables[key_lon][:]   #获取每个像元的经度
    Lat = nc_data_obj.variables[key_lat][:]    #获取每个像元的纬度
    band = np.asarray(nc_data_obj.variables[key_band]).astype(float)  #获取对应波段的像元的值，类型为数组
    print("填充值：",nc_data_obj.variables[key_band])

    #影像的四个角的坐标
    LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] 

    #分辨率计算
    N_Lat = len(Lat)
    if Lon.ndim==1 :
        N_Lon = len(Lon)   #如果Lon为一维的话，就取长度
    else:
        N_Lon = len(Lon[0])   #如果Lon为二维的话，就取宽度
    Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)

    #创建.tif文件
    for i in range(band.shape[0]):
        driver = gdal.GetDriverByName('GTiff')   # 创建驱动
        arr1 = band[i,:,:]                   # 获取不同时间段的数据
        out_tif_name = out_path+os.sep+ data.split('\\')[-1][:4]+ '_'+str(i) +'.tif'
        out_tif = driver.Create(out_tif_name,N_Lon,N_Lat,1,gdal.GDT_Float32)
        # 设置影像的显示范围
        #-Lat_Res一定要是-的
        geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
        out_tif.SetGeoTransform(geotransform)

        #获取地理坐标系统信息，用于选取需要的地理坐标系统
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(coord)                               # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
        out_tif.SetProjection(srs.ExportToWkt())               # 给新建图层赋予投影信息

        #更改异常值
        arr1[arr1[:, :]> 1000000] = -32767

        #数据写出
        if arr1.ndim==2:     #如果本来就是二维数组就不变
            a = arr1[:,:]
        else:                #将三维数组转为二维
            a = arr1[0,:,:]
        out_tif.GetRasterBand(1).WriteArray(a)
        out_tif.GetRasterBand(1).SetNoDataValue(-32767)
        out_tif.FlushCache() # 将数据写入硬盘
        del out_tif       # 注意必须关闭tif文件

'''这个函数部分不需要更改'''

def nc4_to_tif(Input_folder,end_name = 'nc4'):
    Output_folder = os.path.split(Input_folder)[0]  + os.sep+ 'out_' + os.path.split(Input_folder)[-1]
    # 读取所有nc数据
    data_list = glob.glob(Input_folder + os.sep + '*' + end_name)
    print("输入位置为: ",Input_folder)
    print("被读取的nc文件有：",data_list)
    if os.path.exists(Output_folder):
        shutil.rmtree(Output_folder)          #如果文件夹存在就删除
    os.makedirs(Output_folder)            #再重建，这样就不用运行之后又要删了之后再运行了，可以一直运行
    for i in range(len(data_list)):
        dat = data_list[i]
        NC_to_tiffs(data = dat,out_path = Output_folder)
        print (dat + '----转tif成功')
    print(f"输入位置为: {Input_folder}")
    print(f'输出位置为: {Output_folder}')

'''输入路径不能有中文字符----------比如放在D盘中(目前我发现只有有多时间序列的nc或nc4文件会有这个问题，而单时间序列的就不会，这个可以留给大家一起讨论讨论------)'''
nc4_to_tif(Input_folder = r'D:\nc4\nc4',end_name='nc4')   #用户需要输入 ：nc文件所放的文件夹的路径，默认输出至同级目录中，名为'out_...'

在这里，我想补充几点（可能代码的注释里面讲的不是很清楚）：

1.如果想直接使用这个代码的话，只需要修改：

nc4_to_tif(Input_folder = r'D:\nc4\nc4',end_name='nc4')   #用户需要输入 ：nc文件所放的文件夹的路径，默认输出至同级目录中，名为'out_...'

里面的Input_folder的值，这里 r'....' 的意思是防止转义，最好也不要更改。

2.这里用了自动创建文件夹和删除文件夹，这样一来就可以无限次地运行，避免了每运行一次，想再次运行的话，又得重新删除文件夹，用到的代码在这：

if os.path.exists(Output_folder):
        shutil.rmtree(Output_folder)          #如果文件夹存在就删除
    os.makedirs(Output_folder)            #再重建，这样就不用运行之后又要删了之后再运行了，可以一直运行

3.如果大家按照要求运行的话（路径没有中文字符），会出来如下结果：

这里是需要您从 “基本属性” 这里的提示中获取您想要转换为tif数据的波段信息，像这里，我需要的是ndvi这个波段的数据，那我就输入 “ndvi”

点击回车，它就会继续运行，直到输出：

这样就表示输出完成，并且会把输出的路径都给你显示出来，这里我的输出路径为：“D:\nc4\out_nc4”，所以我就可以直接复制，粘贴到搜索目录里面去找这些文件的位置（默认是放在与您输入路径同一级的目录下，名称为 “out_” + “输入的文件名”）。

像：

这里应该都没毛病吧~~~~~~~~

（如果想看代码里面的具体的算法，请看上述的代码的内容以及注释~~~~~~~~）

二、nc_to_tif（多时间序列）

其实这里要说明的话与上面没有什么不同，只是数据由nc4数据变为了nc数据，还有就是代码里面的内容有所不同，操作的话还是一样，一样的~~~可以直接使用，但是如果想深入学习的话，还是得详细看代码哦，里面的注释也是很详细的~~~~~，这里我就不多赘述了~~~~~~，直接上代码

代码如下（示例）：

# -*- coding: utf-8 -*-

import numpy as np
import netCDF4 as nc
from osgeo import gdal,osr,ogr
import glob
import os
from zipfile import ZipFile
import shutil
band_name = ''
def NC_to_tiffs(data,out_path):
    '''
    这个函数里面有些地方还是可能需要更改,像：
    coord(坐标系)
    '''
    coord = 4326            #坐标系，["EPSG","4326"]，默认为4326
    nc_data_obj = nc.Dataset(data)
    #print(nc_data_obj,type(nc_data_obj)) # 了解nc的数据类型，<class 'netCDF4._netCDF4.Dataset'>
    #print(nc_data_obj.variables)      #了解nc数据的基本信息
    key=list(nc_data_obj.variables.keys())            #获取时间，经度，纬度，波段的名称信息，这些可能是不一样的
    print('基础属性为  ',key)
    lon_size = [i for i,x in enumerate(key) if (x.find('lon'.upper())!=-1 or x.find('lon'.lower())!=-1)][0]   #模糊查找属于经度的名称
    lat_size = [i for i,x in enumerate(key) if (x.find('lat'.upper())!=-1 or x.find('lat'.lower())!=-1)][0]   #模糊查找属于纬度的名称
    time_size = [i for i,x in enumerate(key) if (x.find('ime'.upper())!=-1 or x.find('ime'.lower())!=-1)][0]  #模糊查找属于时间的名称
    global band_name
    if band_name == '':
        band_name = input("请输入您想要输出的波段的名字（您可以从'基础属性中得来',不用加上引号）———————:")   #这里是从用户那传入一个波段的字符串，因为nc4的数据比nc要复杂，所以要让用户确定波段的名字
    band_size = [i for i,x in enumerate(key) if (x.find(str(band_name).upper())!=-1 or x.find(str(band_name).lower())!=-1)][0]
    key_band = key[band_size]            #获取波段的名称
    time_name= key[time_size]  #获取时间的名称
    key_lon = key[lon_size]      #获取经度的名称
    key_lat = key[lat_size]      #获取纬度的名称
    Lon = nc_data_obj.variables[key_lon][:]   #获取每个像元的经度
    Lat = nc_data_obj.variables[key_lat][:]    #获取每个像元的纬度
    time = nc_data_obj.variables[time_name]
    times = nc.num2date(time[:],time.units)  # 时间的格式转换,得到一个数组
    band = np.asarray(nc_data_obj.variables[key_band]).astype(float)  #获取对应波段的像元的值，类型为数组
    print("填充值：",nc_data_obj.variables[key_band])

    #影像的四个角的坐标
    LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] 

    #分辨率计算
    N_Lat = len(Lat)
    if Lon.ndim==1 :
        N_Lon = len(Lon)   #获取长度
    else:
        N_Lon = len(Lon[0])
    Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)

    #创建.tif文件
    for i in range(band.shape[0]):
        # strftime() 格式化datetime 对象
        dt = times[i].strftime("%Y-%m")
        driver = gdal.GetDriverByName('GTiff')   # 创建驱动
        arr1 = band[i,:,:]                   # 获取不同时间段的数据
        out_tif_name = out_path+os.sep+ data.split('\\')[-1][:-3]+ dt +'.tif'
        out_tif = driver.Create(out_tif_name,N_Lon,N_Lat,1,gdal.GDT_Float32) 

        # 设置影像的显示范围
        #-Lat_Res一定要是-的
        geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
        out_tif.SetGeoTransform(geotransform)

        #获取地理坐标系统信息，用于选取需要的地理坐标系统
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(coord)                       # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
        out_tif.SetProjection(srs.ExportToWkt())               # 给新建图层赋予投影信息

        #更改异常值
        arr1[arr1[:, :]> 1000000] = -32767

        #数据写出
        if arr1.ndim==2:     #如果本来就是二维数组就不变
            a = arr1[:,:]
        else:     #将三维数组转为二维
            a = arr1[0,:,:]    

        reversed_arr = a[::-1]        #这里是需要倒置一下的     横重要的！！！！！！！！！！！
        out_tif.GetRasterBand(1).WriteArray(reversed_arr)
        out_tif.GetRasterBand(1).SetNoDataValue(-32767)
        out_tif.FlushCache() # 将数据写入硬盘
        del out_tif       # 注意必须关闭tif文件

def nc_to_tif(Input_folder,end_name = 'nc'):
    Output_folder = os.path.split(Input_folder)[0]  + 'out_' + os.path.split(Input_folder)[-1]
    # 读取所有nc数据
    data_list = glob.glob(Input_folder + os.sep + '*' + end_name)
    print("输入位置为: ",Input_folder)
    print("被读取的nc文件有：",data_list)
    #if not os.path.isdir(Output_folder):           #如果输出路径,不存在就创建
    if os.path.exists(Output_folder):
        shutil.rmtree(Output_folder)          #如果文件夹存在就删除
    os.makedirs(Output_folder)            #再重建，这样就不用运行之后又要删了之后再运行了
    for i in range(len(data_list)):
        dat = data_list[i]
        NC_to_tiffs(data = dat,out_path = Output_folder)
        print (dat + '-----转tif成功')
    print(f"输入位置为: {Input_folder}")
    print(f'输出位置为: { Output_folder}')

'''输入路径不能有中文字符----------比如放在D盘中(目前我发现只有有多时间序列的nc或nc4文件才会有这个问题，，单时间序列的nc文件不会出现这样的问题)'''
nc_to_tif(Input_folder = r'D:\spei',end_name='nc')   #用户需要输入 ：nc文件所放的文件夹的路径，默认输出至同一上级目录中

三、nc_to_tif（单时间序列）

这里的要说明的话也和上面一样一样的，所以我就~~~~~~哈哈哈不说太多了哈，不过这里需要用户进行的操作要更少一点。这里是不需要用户传入波段的信息，直接修改下文件的输入路径，就可以直接输出了，并且这里的文件路径可以不用再管是否有中文字符，比较方便哦~~~~~，具体代码的细节都在注释里了哦，爱学习的兄弟可以看看哦~~~~~~~~~~

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

def NC_to_tiffs(data,out_path):
    '''
    这个函数里面有些地方还是可能需要更改
    coord(坐标系)
    '''
    coord = 4326            #坐标系，["EPSG","4326"]，默认为4326
    nc_data_obj = nc.Dataset(data)
    #print(nc_data_obj,type(nc_data_obj)) #了解nc数据的数据类型，<class 'netCDF4._netCDF4.Dataset'>
    #print(nc_data_obj.variables)          #了解nc数据的基本信息
    key=list(nc_data_obj.variables.keys())    #获取时间，经度，纬度，波段的名称信息，这些可能是不一样的
    print('基础属性为',key)
    lon_size = [i for i,x in enumerate(key) if (x.find('lon'.upper())!=-1 or x.find('lon'.lower())!=-1)][0]   #模糊查找属于经度的名称，还在更新.....
    lat_size = [i for i,x in enumerate(key) if (x.find('lat'.upper())!=-1 or x.find('lat'.lower())!=-1)][0]   #模糊查找属于纬度的名称，还在更新.....
    key_band = key[len(key)-1]            #获取波段的名称     目前发现都是放在最后
    key_lon = key[lon_size]      #获取经度的名称
    key_lat = key[lat_size]      #获取纬度的名称
    Lon = nc_data_obj.variables[key_lon][:]#获取每个像元的经度,类型为数组
    Lat = nc_data_obj.variables[key_lat][:]#获取每个像元的纬度，类型为数组
    Band = np.asarray(nc_data_obj.variables[key_band])  #获取对应波段的像元的值，类型为数组
    #影像的四个角的坐标
    LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] 

    #分辨率计算
    N_Lat = len(Lat)
    N_Lon = len(Lon[0])
    Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)

    #创建.tif文件
    driver = gdal.GetDriverByName('GTiff')
    out_tif_name = out_path+os.sep+ data.split('\\')[-1][:-3]+ '.tif'
    out_tif = driver.Create(out_tif_name,N_Lon,N_Lat,1,gdal.GDT_Float32) 

    # 设置影像的显示范围
    #-Lat_Res一定要是-的
    geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
    out_tif.SetGeoTransform(geotransform)

    #获取地理坐标系统信息，用于选取需要的地理坐标系统
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(coord)                       # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
    out_tif.SetProjection(srs.ExportToWkt())               # 给新建图层赋予投影信息

    #更改异常值
    Band[Band[:, :]> 100000] = -32767

    #数据写出
    if Band.ndim==2:     #如果本来就是二维数组就不变
        a = Band[:,:]
    else:       #将三维数组转为二维
        a = Band[0,:,:]
    reversed_arr = a[::-1]      #这里是需要倒置一下的        #这个很重要！！！！！！
    out_tif.GetRasterBand(1).WriteArray(reversed_arr)
    out_tif.GetRasterBand(1).SetNoDataValue(-32767)
    out_tif.FlushCache() # 将数据写入硬盘
    del out_tif       # 注意必须关闭tif文件

def nc_to_tif(Input_folder):
    Output_folder = os.path.split(Input_folder)[0] +os.sep + 'out_' + os.path.split(Input_folder)[1]
    # 读取所有nc数据
    data_list = glob.glob(Input_folder + os.sep + '*.nc')
    print("输入位置为: ",Input_folder)
    print("被读取的nc文件有：",data_list)
#     if not os.path.isdir(Output_folder):
    if os.path.exists(Output_folder):
        shutil.rmtree(Output_folder)          #如果文件夹存在就删除
    os.makedirs(Output_folder)            #再重建，这样就不用运行之后又要删了之后再运行了
    for i in range(len(data_list)):
        dat = data_list[i]
        NC_to_tiffs(data = dat,out_path = Output_folder)
        print (dat + '-----转tif成功')
    print(f"输入位置为: {Input_folder}")
    print("--------------------------")
    print(f'输出位置为: { Output_folder}')

'''#用户需要输入 ：nc文件所放的文件夹的路径，默认输出至同一上级目录中'''  

nc_to_tif(Input_folder = r'D:\nc\real\T2')

总结

还有，还有，还有，这里有几个小坑以及心得我是想我跟大家进行分享de~~~~

1.nc4跟nc的差别在于nc4的数据结构比nc要复杂，内容更丰富，所以转为tif的时候要考虑的东西也更多~~~~~~~~~

2.多时间序列和单时间序列的nc或nc4数据处理成tif形式的方式也不太一样，多时间序列的话要考虑时间因素，单时间是不需要考虑时间因素的，虽然也有时间，但是时间段只有1个，多时间序列的话要根据时间段来进行输出的命名，所以这里也是需要考虑的~~~~~~~~

3.这个是最重要的：就是nc4的数据它是不需要将数据进行颠倒一下的，而nc的数据是需要颠倒的，这个真的是我苦苦发现的，之前也犯了很多很多的错，网上也没有具体的说明，但是这个坑我在代码里面是有说明哦，注释也很详细，所以，如果把上面的代码运行的好的话是不会发生数据颠倒的情况的哦~~~~~~~~~~

到此这篇关于使用python进行nc转tif的3种情况的文章就介绍到这了,更多相关python进行nc转tif内容请搜索我们以前的文章或继续浏览下面的相关文章希望大家以后多多支持我们！