python进行nc转tif

时间:2022-05-23 HeYang.hello 人气:0

前言

本人是位大二在读在校学生，专业为地理信息科学，因跟老师一起做项目，所以有幸接触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的数据是需要颠倒的，这个真的是我苦苦发现的，之前也犯了很多很多的错，网上也没有具体的说明，但是这个坑我在代码里面是有说明哦，注释也很详细，所以，如果把上面的代码运行的好的话是不会发生数据颠倒的情况的哦~~~~~~~~~~

加载全部内容