文章目錄
- 前言
- 一、nc4_to_tif(多時間序列)
- 二、nc_to_tif(多時間序列)
- 三、nc_to_tif(單時間序列)
- 總結
前言
本人是位大二在讀在校學生,專業為地理資訊科學,因跟老師一起做專案,所以有幸接觸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的資料是需要顛倒的,這個真的是我苦苦發現的,之前也犯了很多很多的錯,網上也沒有具體的說明,但是這個坑我在代碼里面是有說明哦,注釋也很詳細,所以,如果把上面的代碼運行的好的話是不會發生資料顛倒的情況的哦~~~~~~~~~~
最后,我想說就是,這是小生我第一次發博客,雖然做專案也有一段時間了,但是也沒發過什么博客,所以,我希望呢~大家能給小生我鼓鼓勵,就比如點點贊,收收藏,評評論是吧,哈哈哈,感謝感謝,大家一起加油哦!!!!~~~~~~~~~~~~~~
轉載請註明出處,本文鏈接:https://www.uj5u.com/houduan/396556.html
標籤:python
上一篇:如何獲得串列串列的列式最小值?
下一篇:使用泛型進行介面向下轉換


