首頁 > 軟體

使用python進行nc轉tif的3種情況解決

2022-03-01 10:00:20

前言

本人是位大二在讀在校學生，專業為地理資訊科學，因跟老師一起做專案，所以有幸接觸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:nc4nc4',end_name='nc4')   #使用者需要輸入 ：nc檔案所放的資料夾的路徑，預設輸出至同級目錄中，名為'out_...'

在這裡，我想補充幾點（可能程式碼的註釋裡面講的不是很清楚）：

1.如果想直接使用這個程式碼的話，只需要修改：

nc4_to_tif(Input_folder = r'D:nc4nc4',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:nc4out_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:ncrealT2')

總結

還有，還有，還有，這裡有幾個小坑以及心得我是想我跟大家進行分享de~~~~

1.nc4跟nc的差別在於nc4的資料結構比nc要複雜，內容更豐富，所以轉為tif的時候要考慮的東西也更多~~~~~~~~~

2.多時間序列和單時間序列的nc或nc4資料處理成tif形式的方式也不太一樣，多時間序列的話要考慮時間因素，單時間是不需要考慮時間因素的，雖然也有時間，但是時間段只有1個，多時間序列的話要根據時間段來進行輸出的命名，所以這裡也是需要考慮的~~~~~~~~

3.這個是最重要的：就是nc4的資料它是不需要將資料進行顛倒一下的，而nc的資料是需要顛倒的，這個真的是我苦苦發現的，之前也犯了很多很多的錯，網上也沒有具體的說明，但是這個坑我在程式碼裡面是有說明哦，註釋也很詳細，所以，如果把上面的程式碼執行的好的話是不會發生資料顛倒的情況的哦~~~~~~~~~~

到此這篇關於使用python進行nc轉tif的3種情況的文章就介紹到這了,更多相關python進行nc轉tif內容請搜尋it145.com以前的文章或繼續瀏覽下面的相關文章希望大家以後多多支援it145.com！


相關文章

  • 德漢姆中學就讀條件 德漢姆中學就讀條件

    <em>Mac</em>Book项目 2009年学校开始实施<em>Mac</em>Book项目,所有师生配备一本<em>Mac</em>Book,并同步更新了校园无线网络。学校每周进行电脑技术更新,每月发送技术支持资料,极大改变了教学及学习方式。因此2011

    2021-06-01 09:32:01

  • 黑科技的輕簡出行，告別「電量焦慮症」：Anker氮化鎵超能充系列 黑科技的輕簡出行，告別「電量焦慮症」：Anker氮化鎵超能充系列

    综合看Anker超能充系列的性价比很高，并且与不仅和iPhone12/苹果<em>Mac</em>Book很配，而且适合多设备充电需求的日常使用或差旅场景，不管是安卓还是Switch同样也能用得上它，希望这次分享能给准备购入充电器的小伙伴们有所

    2021-06-01 09:31:42

  • 吳亦凡廠牌首秀，L4WUDU居然忘詞了，Rapper和明星比還是有差距！ 吳亦凡廠牌首秀，L4WUDU居然忘詞了，Rapper和明星比還是有差距！

    除了L4WUDU与吴亦凡已经多次共事，成为了明面上的厂牌成员，吴亦凡还曾带领20XXCLUB全队参加2020年的一场音乐节，这也是20XXCLUB首次全员合照，王嗣尧Turbo、陈彦希Regi、<em>Mac</em> Ova Seas、林渝植等人全部出场。然而让

    2021-06-01 09:31:34

  • IPFS、Chia、Bzz和ICP挖礦該怎麼選？哪一個更有優勢？ IPFS、Chia、Bzz和ICP挖礦該怎麼選？哪一個更有優勢？

    目前应用IPFS的机构：1 谷歌<em>浏览器</em>支持IPFS分布式协议 2 万维网 (历史档案博物馆)数据库 3 火狐<em>浏览器</em>支持 IPFS分布式协议 4 EOS 等数字货币数据存储 5 美国国会图书馆，历史资料永久保存在 IPFS 6 加

    2021-06-01 09:31:24

  • 有哪些事是買了雪佛蘭後才知道的？美事偷著樂，開拓者車主隨聊 有哪些事是買了雪佛蘭後才知道的？美事偷著樂，開拓者車主隨聊

    开拓者的车机是兼容苹果和<em>安卓</em>，虽然我不怎么用，但确实兼顾了我家人的很多需求：副驾的门板还配有解锁开关，有的时候老婆开车，下车的时候偶尔会忘记解锁，我在副驾驶可以自己开门：第二排设计很好，不仅配置了一个很大的

    2021-06-01 09:30:48

  • iPhone12在618最新定價，跌價幅度超過1400元，還等iPhone13嗎？ iPhone12在618最新定價，跌價幅度超過1400元，還等iPhone13嗎？

    不仅是<em>安卓</em>手机，苹果手机的降价力度也是前所未有了，iPhone12也“跳水价”了，发布价是6799元，如今已经跌至5308元，降价幅度超过1400元，最新定价确认了。iPhone12是苹果首款5G手机，同时也是全球首款5nm芯片的智能机，它

    2021-06-01 09:30:45

IT145.com E-mail:sddin#qq.com