基于Python實(shí)現(xiàn)nc批量轉(zhuǎn)tif格式

更新時(shí)間：2022年08月12日 14:27:24 作者：NCR1008

做項(xiàng)目有時(shí)會(huì)運(yùn)用到netCDF格式的氣象數(shù)據(jù)，而ArcGIS中需要用柵格影像進(jìn)行處理，對(duì)于較多的文件，ArcGIS一個(gè)個(gè)手動(dòng)轉(zhuǎn)換過(guò)于繁瑣，因此我們采用Python進(jìn)行轉(zhuǎn)換，下面就是Python實(shí)現(xiàn)nc批量轉(zhuǎn)tif格式的示例代碼，希望對(duì)你有所幫助

由于做項(xiàng)目需要運(yùn)用到netCDF格式的氣象數(shù)據(jù)，而ArcGIS中需要用柵格影像進(jìn)行處理，對(duì)于較多的文件，ArcGIS一個(gè)個(gè)手動(dòng)轉(zhuǎn)換過(guò)于繁瑣，因此我們采用Python進(jìn)行轉(zhuǎn)換，當(dāng)然也可以采用matlab進(jìn)行轉(zhuǎn)換。

首先需要安裝下面幾個(gè)庫(kù)：

import os
import netCDF4 as nc
import numpy as np
from osgeo import gdal, osr, ogr
import glob

我們可以在下面網(wǎng)址中尋找對(duì)應(yīng)python安裝版本的安裝包，下載后，在pycharm控制臺(tái)中直接安裝即可。例如pip install netCDF4-1.5.8-cp39-cp39-

win_amd64.whl

https://www.lfd.uci.edu/~gohlke/pythonlibs/

安裝之后即可進(jìn)行轉(zhuǎn)換：

def nc2tif(data, Output_folder):
    tmp_data = nc.Dataset(data)  # 利用.Dataset()方法讀取nc數(shù)據(jù)
    print('tmp_data', tmp_data)
 
    Lat_data = tmp_data.variables['lat'][:]
    Lon_data = tmp_data.variables['lon'][:]
    # print(Lat_data)
    # print(Lon_data)
 
    tmp_arr = np.asarray(tmp_data.variables['temp'])
 
    # 影像的左上角&右下角坐標(biāo)
    Lonmin, Latmax, Lonmax, Latmin = [Lon_data.min(), Lat_data.max(), Lon_data.max(), Lat_data.min()]
    # print(Lonmin, Latmax, Lonmax, Latmin)
 
    # 分辨率計(jì)算
    Num_lat = len(Lat_data)  # 5146
    Num_lon = len(Lon_data)  # 7849
    Lat_res = (Latmax - Latmin) / (float(Num_lat) - 1)
    Lon_res = (Lonmax - Lonmin) / (float(Num_lon) - 1)
    # print(Num_lat, Num_lon)
    # print(Lat_res, Lon_res)
 
    for i in range(len(tmp_arr[:])):
        # i=0,1,2,3,4,5,6,7,8,9,...
        # 創(chuàng)建tif文件
        driver = gdal.GetDriverByName('GTiff')
        out_tif_name = Output_folder + '\\' + data.split('\\')[-1].split('.')[0] + '_' + str(i + 1) + '.tif'
        out_tif = driver.Create(out_tif_name, Num_lon, Num_lat, 1, gdal.GDT_Int16)
 
        # 設(shè)置影像的顯示范圍
        # Lat_re前需要添加負(fù)號(hào)
        geotransform = (Lonmin, Lon_res, 0.0, Latmax, 0.0, -Lat_res)
        out_tif.SetGeoTransform(geotransform)
 
        # 定義投影
        prj = osr.SpatialReference()
        prj.ImportFromEPSG(4326)  # WGS84
        out_tif.SetProjection(prj.ExportToWkt())
 
        # 數(shù)據(jù)導(dǎo)出
        out_tif.GetRasterBand(1).WriteArray(tmp_arr[i])  # 將數(shù)據(jù)寫(xiě)入內(nèi)存，此時(shí)沒(méi)有寫(xiě)入到硬盤(pán)
        out_tif.FlushCache()  # 將數(shù)據(jù)寫(xiě)入到硬盤(pán)
        out_tif = None  # 關(guān)閉tif文件
 
def main():
    Input_folder = r"E:\competition\航天宏圖\2-meter air temperature_CMFD\Data_forcing_01yr_010deg\\"
    # Input_folder = r"E:\competition\航天宏圖\2-meter air temperature_CMFD\Data_forcing_01mo_010deg\\"
    Output_folder = r"E:\competition\航天宏圖\2-meter air temperature_CMFD\tif\\"
 
    # 讀取所有數(shù)據(jù)
    data_list = glob.glob(os.path.join(Input_folder + '*.nc'))
    print(data_list)
 
    for i in range(len(data_list)):
        data = data_list[i]
        nc2tif(data, Output_folder)
        print(data + '轉(zhuǎn)tif成功')

值得注意的是，tmp_arr = np.asarray(tmp_data.variables['temp'])中的temp可以根據(jù)需要轉(zhuǎn)換的波段進(jìn)行選擇，我們可以在讀取數(shù)據(jù)之后print一下，找到對(duì)應(yīng)的波段進(jìn)行替換即可。如下圖中我們應(yīng)該選擇的就是temp。

完成上述步驟即可得到所需的tif圖像：

在上述代碼中，經(jīng)過(guò)處理的影像是倒置的，可能是處理過(guò)程中仿射矩陣讀寫(xiě)錯(cuò)誤導(dǎo)致的。因此我們可以在寫(xiě)入影像的時(shí)候，進(jìn)行影像的垂直鏡像操作即可：WriteArray(ndvi_arr_float[i][::-1])

def NC_to_tiffs(data, Output_folder):
    nc_data_obj = nc.Dataset(data)
    Lon = nc_data_obj.variables['lon'][:]
    Lat = nc_data_obj.variables['lat'][:]
    ndvi_arr = np.asarray(nc_data_obj.variables['temp'])  
    ndvi_arr_float = ndvi_arr.astype(float) / 10000  之間
    # 影像的左上角和右下角坐標(biāo)
    LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]
    # 分辨率計(jì)算
    N_Lat = len(Lat)
    N_Lon = len(Lon)
    Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1)
    for i in range(len(ndvi_arr[:])):
        driver = gdal.GetDriverByName('GTiff')
        out_tif_name = Output_folder + '\\' + data.split('\\')[-1].split('.')[0] + '_' + str(i + 1) + '.tif'
        out_tif = driver.Create(out_tif_name, N_Lon, N_Lat, 1, gdal.GDT_Float32)
       
        geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
        out_tif.SetGeoTransform(geotransform)
      
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326) 
        out_tif.SetProjection(srs.ExportToWkt()) 
        # 數(shù)據(jù)寫(xiě)出
        out_tif.GetRasterBand(1).WriteArray(ndvi_arr_float[i][::-1])  # 將數(shù)據(jù)寫(xiě)入內(nèi)存，此時(shí)沒(méi)有寫(xiě)入硬盤(pán) 此處[::-1]用于圖像的垂直鏡像對(duì)稱，避免圖像顛倒
        out_tif.FlushCache()  # 將數(shù)據(jù)寫(xiě)入硬盤(pán)
        out_tif = None  # 注意必須關(guān)閉tif文件

這樣便可以得到正確輸出的影像：

當(dāng)然，我們除了在寫(xiě)入時(shí)做垂直鏡像操作之外，還可以利用python對(duì)圖像進(jìn)行幾何變換來(lái)實(shí)現(xiàn)翻轉(zhuǎn)。具體代碼如下：

圖像水平翻轉(zhuǎn)：

#  圖像水平翻轉(zhuǎn)
    im_data_hor = np.flip(im_data, axis=2)
    hor_path = train_image_path + "\\" + str(tran_num) + imageList[i][-4:]
    writeTiff(im_data_hor, im_geotrans, im_proj, hor_path)

標(biāo)簽水平翻轉(zhuǎn)：

 #  標(biāo)簽水平翻轉(zhuǎn)
    Hor = cv2.flip(label, 1)
    hor_path = train_label_path + "\\" + str(tran_num) + labelList[i][-4:]
    cv2.imwrite(hor_path, Hor)
    tran_num += 1

圖像垂直翻轉(zhuǎn)：

 #  圖像垂直翻轉(zhuǎn)
    im_data_vec = np.flip(im_data, axis=1)
    vec_path = train_image_path + "\\" + str(tran_num) + imageList[i][-4:]
    writeTiff(im_data_vec, im_geotrans, im_proj, vec_path)

標(biāo)簽垂直翻轉(zhuǎn)：

 #  標(biāo)簽垂直翻轉(zhuǎn)
    Vec = cv2.flip(label, 0)
    vec_path = train_label_path + "\\" + str(tran_num) + labelList[i][-4:]
    cv2.imwrite(vec_path, Vec)
    tran_num += 1

圖像鏡像翻轉(zhuǎn)：

 #  圖像對(duì)角鏡像
    im_data_dia = np.flip(im_data_vec, axis=2)
    dia_path = train_image_path + "\\" + str(tran_num) + imageList[i][-4:]
    writeTiff(im_data_dia, im_geotrans, im_proj, dia_path)

標(biāo)簽鏡像翻轉(zhuǎn)：

 #  標(biāo)簽對(duì)角鏡像
    Dia = cv2.flip(label, -1)
    dia_path = train_label_path + "\\" + str(tran_num) + labelList[i][-4:]
    cv2.imwrite(dia_path, Dia)
    tran_num += 1

若是輸出路徑的文件夾沒(méi)有建立好，則會(huì)報(bào)如下錯(cuò)誤。當(dāng)然，為了減少工作量，也可以定義一個(gè)函數(shù)，如果路徑不存在則自動(dòng)創(chuàng)建，就可以解決這個(gè)問(wèn)題。

到此這篇關(guān)于基于Python實(shí)現(xiàn)nc批量轉(zhuǎn)tif格式的文章就介紹到這了,更多相關(guān)Python nc轉(zhuǎn)tif內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: