区域生态质量评价办法（试行）——原理及算法实现（Python）

目录

适用范围

指标体系

1.生态格局

（1）生态用地面积比指数

（2）海洋自然岸线保有指数

（3）生境质量指数

适用范围

本办法规定了区域生态质量评价的指标体系，数据要求和评价办法。适用于县级及以上区域生态质量现状和趋势的综合评价。

指标体系

指标体系包含了生态格局、生态功能、生物多样性和生态胁迫四个一级指标，下设11个二级指标、18个三级指标。

一级指标	二级指标	三级指标	备注
生态格局	生态组分	生态用地面积比指数
		海洋自然岸线保有指数	沿海县域
	生态结构	生态保护红线面积比指数
		生境质量指数
		重要生态空间连通度指数
生态功能	水土保持	水土保持指数	水土保持类型国家重点生态功能区县域
	水源涵养	水源涵养指数	水源涵养类型国家重点生态功能区县域
	防风固沙	防风固沙指数	防风固沙类型国家重点生态功能区县域
	生态宜居	建成区绿地率指数	地级及以上城市建成区
		建成区公园绿地可达指数
	生态活力	植被覆盖指数	其他县域
		水网密度指数
生物多样性	生物保护	重点保护生物指数
	重要生物功能群	只是生物类群生命力指数
		原生功能群种占比指数
生态胁迫	人为胁迫	陆域开发干扰指数
		陆域开发强度指数	沿海县域
	自然胁迫	自然灾害受灾指数

1.生态格局

（1）生态用地面积比指数

使用数据：以二调或三调矢量数据为基础，根据生态用地的类型，筛选对应生态用地计算面积，并利用权重参数计算最终的生态用地面积比指数。

def EL(data):
    
    forest = np.where((data[:,2]<24))
    forest = np.array(forest,dtype=int)
    data_forest = data[forest[0][:],:]
    area_f = cal_area(data_forest)

    grass = np.where((data[:,2]<34) & (data[:,2]>30))
    grass = np.array(grass,dtype=int)
    data_grass = data[grass[0][:],:]
    area_g = cal_area(data_grass)

    water = np.where((data[:,2]>40) & (data[:,2]<47))
    water = np.array(water,dtype=int)
    data_water = data[water[0][:],:]
    area_w = cal_area(data_water)

    sand = np.where((data[:,2]==61))
    sand = np.array(sand,dtype=int)
    data_sand = data[sand[0][:],:]
    area_sd = cal_area(data_sand)

    o_forest = np.where((data[:,2]==24))
    o_forest = np.array(o_forest,dtype=int)
    data_o_forest = data[o_forest[0][:],:]
    area_of = cal_area(data_o_forest)

    water_sk = np.where((data[:,2]==43))
    water_sk = np.array(water_sk,dtype=int)
    data_water_sk = data[water_sk[0][:],:]
    area_sk = cal_area(data_water_sk)

    farm_st = np.where((data[:,2]<115) & (data[:,2]>110))
    farm_st = np.array(farm_st,dtype=int)
    data_farm_st = data[farm_st[0][:],:]
    area_st = cal_area(data_farm_st) 

    farm_hd = np.where((data[:,2]<125) & (data[:,2]>120))
    farm_hd = np.array(farm_hd,dtype=int)
    data_farm_hd = data[farm_hd[0][:],:]
    area_hd = cal_area(data_farm_hd)    

    area = cal_area(data)

    result_el = weight['ael']*(area_f+area_g+area_w+area_sd+area_of*weight['other_forest']+area_sk*weight['reservoir']+area_st*weight['paddy_filed']+area_hd*weight['land'])/area
    
return result_el

#通过筛选对应生态用地类型（矢量），并计算面积后套用公式计算生态用地面积比指数

（2）海洋自然岸线保有指数

用于沿海城市海岸线中的自然岸线计算。其中海岸线数据可根据高分遥感影像目视解译获得或从对应部门获取。

#encoding = utf-8
#time: 2022.7.19 11:37:48
#author: Shi-hanning 


import numpy as np 
import arcpy 
import os,sys 

fc = sys.argv[1]
loca = sys.argv[2]

loca = eval(loca)

field = ['xdm','ld2000','length','自然岸线']
values = [row for row in arcpy.da.SearchCursor(fc,field)]
values = np.array(values,dtype='object')

def cal_length(data):
    length = sum(data[:,2])
    return length 

def cal_nature(data):
    ind = np.where(data[:,3]==1)
    ind = np.array(ind,dtype=int)
    data1 = data[ind[0][:],:]
    length = cal_length(data1)
    return length 

result = np.zeros((len(loca),3))
for i in range(len(loca)):
    place = int(loca[i])

    xdm = np.array(values[:,0],dtype='int')
    ind = np.where(xdm==int(place))
    ind = np.array(ind,dtype=int)
    data = values[ind[0][:],:]
    data = np.array(data)

    nature = cal_nature(data)
    length =cal_length(data)
    nonc = 100*nature/length

    result[i][0] = place
    result[i][1] = nature
    result[i][2] = nonc

np.set_printoptions()
out = result.tolist()
print(out)

（3）生境质量指数

评价区域内生态环境质量状况。

field = ['XZDM','Shape_Area']
        values = [row for row in(arcpy.da.SearchCursor(fc,field))]
        values = np.array(values)
        result = np.zeros((len(para),2))
        idx = np.where(values[:,0]==int(place))
        idx = np.array(idx,dtype=int)
        data = values[idx[0][:],:]
        total_area = cal_area(data)
        hqi = weights['Abio']*(weights['forest']*parameter['SF']+
                                weights['grass']*parameter['SG']+
                                weights['water']*parameter['SW']+
                                weights['farmland']*parameter['SC']+
                                weights['urban']*parameter['SB']+
                                weights['other']*parameter['SU'])/total_area

        result[i][0] = parameter['xzdm']
        result[i][1] = hqi



def sb(data):
    city1 = np.where((data[:,2]==51))
    city2 = np.where((data[:,2]==52))
    city3 = np.where((data[:,2]==53))

    city1 = np.array(city1,dtype=int)
    city2 = np.array(city2,dtype=int)
    city3 = np.array(city3,dtype=int)

    data_city1 = data[city1[0][:],:]
    data_city2 = data[city2[0][:],:]
    data_city3 = data[city3[0][:],:]

    area_city1 = cal_area(data_city1)
    area_city2 = cal_area(data_city2)
    area_city3 = cal_area(data_city3)

    #sb = area_city1*0.30 + area_city2*0.40 + area_city3*0.30
    sb = area_city1*weight['u1']+ area_city2*weight['u2'] + area_city3*weight['u3']

    return sb

##SF、GS、SW、SC、SU计算函数同SB，仅需将读取的字段改为对应用地类型的代码即可