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

前言

本人是位大二在读在校学生，专业为地理信息科学，因跟老师一起做项目，所以有幸接触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.如果大家按照要求运行的话（路径没有中文字符），会出来如下结果：

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

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

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

像：

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

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

二、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')

来源：https://blog.csdn.net/weixin_56872334/article/details/122199373