使用python进行nc转tif的3种情况解决
时间：2022-08-27 09:39:38 | 栏目：Python代码 | 点击：次
前言
本人是位大二在读在校学生，专业为地理信息科学，因跟老师一起做项目，所以有幸接触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的数据是需要颠倒的，这个真的是我苦苦发现的，之前也犯了很多很多的错，网上也没有具体的说明，但是这个坑我在代码里面是有说明哦，注释也很详细，所以，如果把上面的代码运行的好的话是不会发生数据颠倒的情况的哦~~~~~~~~~~