python进行nc转tif

软件发布|下载排行|最新软件

当前位置:首页IT学院IT技术

python进行nc转tif

HeYang.hello   2022-05-23 我要评论

前言

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

Copyright 2022 版权所有 软件发布 访问手机版

声明:所有软件和文章来自软件开发商或者作者 如有异议 请与本站联系 联系我们