使用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.如果大傢按照要求運行的話(路徑沒有中文字符),會出來如下結果:

 這裡是需要您從 “基本屬性” 這裡的提示中獲取您想要轉換為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內容請搜索WalkonNet以前的文章或繼續瀏覽下面的相關文章希望大傢以後多多支持WalkonNet!

推薦閱讀: