任務描述

獲取研究區HWSD數據集中：沙含量、淤泥含量、黏土含量、有機碳含量等。

技術流程

1.數據下載
2.影像連接數據庫
3.提取所需土壤圖層

1.數據下載

所需數據在寒區旱區科學數據中心下載HWSD中國土壤數據集
下載的數據包包括：數據說明pdf、img格式、mdb數據庫

2.影像連接數據庫

*比較關鍵的處理步驟為：img+link+dbf數據庫。個人理解：img中只是展示各類土壤的地理分布，像元值代表其土壤分類，用唯一ID值與dbf數據庫的各類土壤屬性相鏈接。具體操作可以參考鏈接

3.提取所需土壤圖層

代碼思路是：逐像元讀取土壤ID→憑借ID讀取excel表中對應土壤屬性值→將屬性值賦值給像元生成新的土壤圖層影像
*土壤屬性表格（自己整理出來）

具體代碼：

# 讀寫遙感數據import gdal import numpy as np from collections import Counterimport xlrd #導入操作excel的庫 import xlwt import osexcel = xlrd.open_workbook('F:\\土壤數據\\白銀WHSD表格數據大范圍.xlsx') table = excel.sheets() [0]#選擇sheet nrows = table.nrows #獲取行數，用來指導循環多少次 print (nrows)filename = (r'F:\\土壤數據\\HWSD_Baiyin_Clip.tif') #文件路徑data = gdal.Open(filename) #打開文件 cols = data.RasterXSize #讀取圖像列數 rows = data.RasterYSize #讀取圖像行數im_geotrans = data.GetGeoTransform() # 仿射矩陣，左上角像素的大地坐標和像素分辨率。 # 共有六個參數，分表代表左上角x坐標；東西方向上圖像的分辨率；如果北邊朝上，地圖的旋轉角度，0表示圖像的行與x軸平行；左上角y坐標； # 如果北邊朝上，地圖的旋轉角度，0表示圖像的列與y軸平行；南北方向上地圖的分辨率。 im_proj = data.GetProjection() # 地圖投影信息 im_data = data.ReadAsArray(0, 0, cols, rows) # 此處讀取整張圖像 # ReadAsArray(<xoff>, <yoff>, <xsize>, <ysize>) # 讀出從(xoff,yoff)開始，大小為(xsize,ysize)的矩陣。 #print (im_data)del data # 釋放內存，如果不釋放，在arcgis，envi中打開該圖像時會顯示文件被占用 print ( cols, rows)list = [] MU_GLOBAL_List = [] T_GRAVEL_List = [] #碎石體積百分比 T_SAND_List = [] #沙含量 T_SILT_List = [] #淤泥含量 T_CLAY_List = [] #粘土含量 T_USDA_TEX_CLASS_List = [] #USDA土壤質地分類 T_REF_BULK_DENSITY_List = [] #USDA土壤質地分類 T_OC_List = [] #有機碳含量 T_PH_H2O_List = [] #酸堿度 T_CEC_CLAY_List = [] #粘性層土壤的陽離子交換能力 T_CEC_SOIL_List = [] #土壤的陽離子交換能力 T_BS_List = [] #基本飽和度 T_TEB_List = [] #交換性鹽基 T_CACO3_List = [] #碳酸鹽或石灰含量 T_CASO4_List = [] #硫酸鹽含量 T_ESP_List = [] #可交換鈉鹽 T_ECE_List = [] #電導率#讀取土壤ID for k in range(1,nrows):MU_GLOBAL = int(table.cell(k, 0).value)MU_GLOBAL_List.append(MU_GLOBAL)T_GRAVEL = int(table.cell(k, 1).value)T_GRAVEL_List.append(T_GRAVEL)T_SAND = int(table.cell(k, 2).value)T_SAND_List.append(T_SAND)T_SILT = int(table.cell(k, 3).value)T_SILT_List.append(T_SILT)T_CLAY = int(table.cell(k, 4).value)T_CLAY_List.append(T_CLAY)T_USDA_TEX_CLASS = int(table.cell(k, 5).value)T_USDA_TEX_CLASS_List.append(T_USDA_TEX_CLASS)T_REF_BULK_DENSITY = float(table.cell(k, 6).value)T_REF_BULK_DENSITY_List.append(T_REF_BULK_DENSITY)T_OC = float(table.cell(k, 7).value)T_OC_List.append(T_OC)T_PH_H2O = float(table.cell(k, 8).value)T_PH_H2O_List.append(T_PH_H2O)T_CEC_CLAY = int(table.cell(k, 9).value)T_CEC_CLAY_List.append(T_CEC_CLAY)T_CEC_SOIL = int(table.cell(k, 10).value)T_CEC_SOIL_List.append(T_CEC_SOIL)T_BS = int(table.cell(k, 11).value)T_BS_List.append(T_BS)T_TEB = float(table.cell(k, 12).value)T_TEB_List.append(T_TEB)T_CACO3 = float(table.cell(k, 13).value)T_CACO3_List.append(T_CACO3)T_CASO4 = float(table.cell(k, 14).value)T_CASO4_List.append(T_CASO4)T_ESP = int(table.cell(k, 15).value)T_ESP_List.append(T_ESP)T_ECE = float(table.cell(k, 16).value)T_ECE_List.append(T_ECE) #print (MU_GLOBAL_List,T_GRAVEL_List,T_SAND_List,T_SILT_List,T_CLAY_List,T_USDA_TEX_CLASS_List,T_REF_BULK_DENSITY_List,T_OC_List, # T_PH_H2O_List,T_CEC_CLAY_List,T_CEC_SOIL_List,T_BS_List,T_TEB_List,T_CACO3_List,T_CASO4_List,T_ESP_List,T_ECE_List)for col in range (rows):for row in range (cols):value = im_data[col][row]list.append(value)for k in range(1, nrows):MU_GLOBAL = int(table.cell(k, 0).value)if value == MU_GLOBAL:NewValue = int(table.cell(k, 2).value) #沙含量#im_data = im_data.astype(np.float) #強行轉換為float格式im_data[col][row] = NewValueelse:continue #print (im_data) #print (list) #result = Counter(list)# 統計不同土壤類型出現的次數，即區域內有多少類型的土壤 #print (result)driver = gdal.GetDriverByName('GTiff') # 數據類型必須有，因為要計算需要多大內存空間 #data = driver.Create('F:\\土壤數據\\有機碳含量.tif', cols, rows, 1,gdal.GDT_Float32 ) #用來輸出浮點型影像 data = driver.Create('F:\\土壤數據\\沙含量.tif', cols, rows) #用來輸出整型影像 data.SetGeoTransform(im_geotrans) # 寫入仿射變換參數 data.SetProjection(im_proj) # 寫入投影 data.GetRasterBand(1).WriteArray(im_data) # 寫入數組數據del data

注：
*im_data = im_data.astype(np.float) #強行轉換為float格式
此句為強行轉化為float格式，因為im_data最開始讀取的img格式為整型，當土壤屬性值為整數時，可以正常讀寫；但當土壤屬性值為浮點型時，則輸出圖像還未整型，值錯誤，所以要進行強行轉化。

*data = driver.Create(‘F:\土壤數據\有機碳含量.tif’, cols, rows, 1,gdal.GDT_Float32 ) #用來輸出浮點型影像
此句與上句配套，以便合適輸出浮點型影像。

總結

以上是生活随笔為你收集整理的土壤HWSD处理流程的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。