土壤HWSD处理流程

任务描述

获取研究区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 Counter

import xlrd #导入操作excel的库
import xlwt
import os

excel = 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] = NewValue
            else:
                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 ) #用来输出浮点型影像
此句与上句配套，以便合适输出浮点型影像。