一.哨兵2號獲取/處理/提取數(shù)據(jù)

1.1 影像處理與下載

這里采用哨兵12號影像估算森林生物量
在GEE上處理和下載2017年的S2L1C級產(chǎn)品，因?yàn)镾2L2A級產(chǎn)品（經(jīng)過大氣校正）量少，沒有2017年的可用產(chǎn)品。

這里需要對S2L1C產(chǎn)品進(jìn)行大氣較正，采用6S模型，運(yùn)行這個庫需要安裝。
可以看這篇文章完成

#導(dǎo)入必要的庫文件
import ee
from Py6S import *
import datetime
import math
import os
import sys
sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin'))
from atmospheric import Atmospheric
ee.Initialize()
import geemap
wanzhou = geemap.geojson_to_ee("./input/wanzhou.json")
roi = wanzhou.geometry().bounds()
geom = wanzhou.geometry().bounds()

采用云概率影像去云

s2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED")
s2_cloud = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY")
point = roidef rmCloudByProbability(image, thread):prob = image.select("probability")return image.updateMask(prob.lte(thread))def scaleImage(image):time_start = image.get("system:time_start")image = image.divide(10000)image = image.set("system:time_start", time_start)return imagedef getMergeImages(primary, secondary):join = ee.Join.inner()filter = ee.Filter.equals(leftField='system:index', rightField='system:index')joinCol = join.apply(primary, secondary, filter)joinCol = joinCol.map(lambda image: ee.Image(image.get("primary")).addBands(ee.Image(image.get("secondary"))))return ee.ImageCollection(joinCol)startDate = "2016-06-01"
endDate = "2016-10-31"
ee.Date(startDate)
ee.Date(endDate)
s2_raw = s2.filterDate(startDate, endDate).filterBounds(point)
s2Imgs1 = s2.filterDate(startDate, endDate).filterBounds(point).map(scaleImage)
s2Imgs2 = s2_cloud.filterDate(startDate, endDate).filterBounds(point)
s2Imgs = getMergeImages(s2Imgs1, s2Imgs2)
s2Imgs = s2Imgs.map(lambda image: rmCloudByProbability(image,40))
s2Img = s2Imgs.median()composite = s2Img.clip(point).toFloat()
rgbVis = {'min': 0.0,'max': 0.3,'bands': ['B4', 'B3', 'B2'],
}
styling = {'color': 'red','fillColor': '00000000'
}
#隨后創(chuàng)建一個 Map 實(shí)例，將柵格和矢量添加到圖層中。
Map = geemap.Map()
Map.addLayer(composite,rgbVis, "S2_2A_wanzhou")
Map.addLayer(wanzhou.style(**styling), {}, '萬州邊界')
Map.centerObject(wanzhou, 9)
Map

toa = composite

采用6S模型對1C級產(chǎn)品進(jìn)行大氣校正

date = ee.Date('2016-07-01')geom = ee.Geometry.Point(107.7632, 30.8239)
S2 = ee.Image(ee.ImageCollection('COPERNICUS/S2').filterBounds(geom).filterDate(date,date.advance(3,'month')).sort('system:time_start').first())

info = S2.getInfo()['properties']
scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)# i.e. Python uses seconds, EE uses milliseconds
solar_z = info['MEAN_SOLAR_ZENITH_ANGLE']
date = ee.Date('2016-07-01')h2o = Atmospheric.water(geom,date).getInfo()
o3 = Atmospheric.ozone(geom,date).getInfo()
aot = Atmospheric.aerosol(geom,date).getInfo()SRTM = ee.Image('CGIAR/SRTM90_V4')# Shuttle Radar Topography mission covers *most* of the Earth
alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom.centroid()).get('elevation').getInfo()
km = alt/1000 # i.e. Py6S uses units of kilometers# Instantiate
s = SixS()# Atmospheric constituents
s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3)
s.aero_profile = AeroProfile.Continental
s.aot550 = aot# Earth-Sun-satellite geometry
s.geometry = Geometry.User()
s.geometry.view_z = 0               # always NADIR (I think..)
s.geometry.solar_z = solar_z        # solar zenith angle
s.geometry.month = scene_date.month # month and day used for Earth-Sun distance
s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance
s.altitudes.set_sensor_satellite_level()
s.altitudes.set_target_custom_altitude(km)def spectralResponseFunction(bandname):"""Extract spectral response function for given band name"""bandSelect = {'B1':PredefinedWavelengths.S2A_MSI_01,'B2':PredefinedWavelengths.S2A_MSI_02,'B3':PredefinedWavelengths.S2A_MSI_03,'B4':PredefinedWavelengths.S2A_MSI_04,'B5':PredefinedWavelengths.S2A_MSI_05,'B6':PredefinedWavelengths.S2A_MSI_06,'B7':PredefinedWavelengths.S2A_MSI_07,'B8':PredefinedWavelengths.S2A_MSI_08,'B8A':PredefinedWavelengths.S2A_MSI_8A,'B9':PredefinedWavelengths.S2A_MSI_09,'B10':PredefinedWavelengths.S2A_MSI_10,'B11':PredefinedWavelengths.S2A_MSI_11,'B12':PredefinedWavelengths.S2A_MSI_12,}return Wavelength(bandSelect[bandname])def toa_to_rad(bandname):"""Converts top of atmosphere reflectance to at-sensor radiance"""# solar exoatmospheric spectral irradianceESUN = info['SOLAR_IRRADIANCE_'+bandname]solar_angle_correction = math.cos(math.radians(solar_z))# Earth-Sun distance (from day of year)doy = scene_date.timetuple().tm_ydayd = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year# conversion factormultiplier = ESUN*solar_angle_correction/(math.pi*d**2)# at-sensor radiancerad = toa.select(bandname).multiply(multiplier)return raddef surface_reflectance(bandname):"""Calculate surface reflectance from at-sensor radiance given waveband name"""# run 6S for this wavebands.wavelength = spectralResponseFunction(bandname)s.run()# extract 6S outputsEdir = s.outputs.direct_solar_irradiance             #direct solar irradianceEdif = s.outputs.diffuse_solar_irradiance            #diffuse solar irradianceLp   = s.outputs.atmospheric_intrinsic_radiance      #path radianceabsorb  = s.outputs.trans['global_gas'].upward       #absorption transmissivityscatter = s.outputs.trans['total_scattering'].upward #scattering transmissivitytau2 = absorb*scatter                                #total transmissivity# radiance to surface reflectancerad = toa_to_rad(bandname)ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif))return ref# # surface reflectance rgb
# b = surface_reflectance('B2')
# g = surface_reflectance('B3')
# r = surface_reflectance('B4')
# ref = r.addBands(g).addBands(b)# all wavebands
output = S2.select('QA60')
for band in ['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10','B11','B12']:
#     print(band)output = output.addBands(surface_reflectance(band))

#隨后創(chuàng)建一個 Map 實(shí)例，將柵格和矢量添加到圖層中。
Map = geemap.Map()
Map.addLayer(output, rgbVis, '萬州')
Map.addLayer(wanzhou.style(**styling), {}, '萬州邊界')
Map.centerObject(wanzhou, 9)
Map

geemap下載數(shù)據(jù)到本地

這里需要安裝geemap庫，具體可以參考這篇文章

#定義下載函數(shù)
def download(image):band_name = image.bandNames().getInfo()band_name_str = str(band_name).replace('[', '').replace(']', '').replace("'", '')work_dir = os.path.join(os.path.expanduser("E:\jupyter\geeDownloads"), 'tif')if not os.path.exists(work_dir):os.makedirs(work_dir)out_tif = os.path.join(work_dir, str(band_name_str)+"_S2_2016_wanzhou.tif")geemap.download_ee_image(image=image,filename=out_tif,region=wanzhou.geometry(),crs="EPSG:4326",scale=10,)return image

#下載原始數(shù)據(jù)download(output)

NDVI

# 選擇 Sentinel-2 的紅波段（B4）和近紅外波段（B8）
red = output.select('B4')
nir = output.select('B8')# 計(jì)算歸一化植被指數(shù)（NDVI）
ndvi = nir.subtract(red).divide(nir.add(red)).rename('ndvi')# 設(shè)置可視化參數(shù)
vis_params = {'min': -1,'max': 1,'palette': ['blue', 'white', 'green']
}
ndvi = ndvi.clip(wanzhou)
# 在地圖上添加 NDVI 圖層
Map.addLayer(ndvi, vis_params, "NDVI_S2_2016")
Map
download(ndvi)

1.2 各種參數(shù)計(jì)算（生物物理變量、植被指數(shù)等）

LAI：葉面積指數(shù)

這幾個參數(shù)都用用snap軟件打開，但是從GEE上下載的數(shù)據(jù)是tif格式的，缺失了頭文件，不能用SNAP處理，我也不知道有什么好的解決方法嗚嗚嗚

參考鏈接

參考文獻(xiàn)：Boegh E,Soegaard H,Broge N,et al.Airborne multispectral data for
quantifying leaf area index,nitrogen concentration,and photosynthetic
efficiency in agriculture[J].Remote Sensing of
Environment,2002,81(2):179-193.

公式：LAI=3.618*EVI-0.118

EVI = output.expression('2.5*((NIR - RED)/(NIR + 6*RED-7.5*BLUE+1))',{'NIR':output.select('B8'),'RED':output.select('B4'),'BLUE':output.select('B2')}
).float().rename('EVI')LAI = output.expression('3.618*EVI-0.118',{'EVI':EVI.select('EVI'),}
).rename('LAI')
LAI=LAI.clip(wanzhou)download(LAI)

FAPAR：吸收的光合有效輻射的分?jǐn)?shù)

參考鏈接

LAI和FAPAR之間存在一定的數(shù)學(xué)關(guān)系，通常使用經(jīng)驗(yàn)公式來描述它們之間的關(guān)系。其中最常用的公式是：

FAPAR = exp（-k * LAI）

其中，k是一個常數(shù)，通常取值在0.5左右。這個公式表明，FAPAR與LAI成指數(shù)反比關(guān)系，即LAI越高，FAPAR越低。

這個公式的基本原理是，LAI越高，表示單位面積內(nèi)植物葉面積越大，可以吸收更多的光能，導(dǎo)致FAPAR值降低。而當(dāng)賴

需要注意的是，這個公式是一種經(jīng)驗(yàn)公式，適用于各種植被類型和環(huán)境條件。在實(shí)際應(yīng)用中，還需要考慮到植被的結(jié)構(gòu)、生長狀態(tài)、光譜特征等因素，以獲得更準(zhǔn)確的LAI和FAPAR估算結(jié)果。

FAPAR = output.expression('exp(-0.51*LAI)',{'LAI':LAI.select('LAI')}
).float().rename('FAPAR')download(FAPAR)

FVC：植被覆蓋率

參考鏈接1

參考鏈接2

FVC = (NDVI -NDVIsoil)/ ( NDVIveg - NDVIsoil)

在ndvi影像值統(tǒng)計(jì)結(jié)果中，最后一列表示對應(yīng)NDVI值的累積概率分布。我們分別取累積概率為5%和90%的NDVI值作為NDVImin和NDVImax

#計(jì)算NDVImin和NDVImax
from osgeo import gdal
import numpy as npinputRaster = gdal.Open(r"E:\jupyter\geeDownloads\NDVI_2016_wanzhou.tif")  # 輸入的遙感影像數(shù)據(jù)
# inputRaster = gdal.Open(r"E:\jupyter\geeDownloads\tif\NDVI.tif")  # 輸入的遙感影像數(shù)據(jù)
ndviBand = inputRaster.GetRasterBand(1)
ndviValues = ndviBand.ReadAsArray()  # 將 NDVI 柵格數(shù)據(jù)轉(zhuǎn)換為 numpy 數(shù)組# 統(tǒng)計(jì) NDVI 值的分布
ndviHist, ndviBins = np.histogram(ndviValues, bins=1000, range=(-1, 1))
ndviCumHist = np.cumsum(ndviHist) / ndviValues.size# 找到累積概率值分別為 0.05 和 0.9 時對應(yīng)的 NDVI 值
ndviBinSize = ndviBins[1] - ndviBins[0]
ndviSoilIndex = np.where(ndviCumHist <= 0.05)[0][-1]
ndviSoil = ndviBins[ndviSoilIndex] + ndviBinSize / 2
ndviVegIndex = np.where(ndviCumHist <= 0.9)[0][-1]
ndviVeg = ndviBins[ndviVegIndex] + ndviBinSize / 2# 輸出統(tǒng)計(jì)結(jié)果
print("ndviSoil: {:.3f}".format(ndviSoil))
print("ndviVeg: {:.3f}".format(ndviVeg))

然后在ArcGis中柵格計(jì)算器中計(jì)算：
Con((“NDVI.tif” >= 0.503) & (“NDVI.tif” <= 0.999), ((“NDVI.tif” - 0.503) / (0.999 - 0.503)), 0)

GEE計(jì)算植被指數(shù)

這是可選項(xiàng)，如果你wifi質(zhì)量好，可以這樣下載，當(dāng)然還是建議本地計(jì)算柵格

1.比值植被指數(shù)RVI =B8/B4

2.紅外植被指數(shù)IPVI=B8/(B8+B4)

3.垂直植被指數(shù)PVI=SIN(45)*B8-COS(45)*B4

4.反紅邊葉綠素指數(shù)IRECI=(B8-B4)/(B8+B4)

5.土壤調(diào)節(jié)植被指數(shù)SAVI=((B8-B4)/(B8+B4+L))(1+L) L取0.5

6.大氣阻抗植被指數(shù)ARVI=B8-(2*B4-B2)/B8+(2B4-B2)

7.特定色素簡單比植被指數(shù)PSSRa=B7/B4

8.Meris陸地葉綠素指數(shù)MTCI=(B6-B5)/(B5-B4)

9.修正型葉綠素吸收比值指數(shù)MCARI=[B5-B4]-0.2*(B5-B3)]*(B5-B4)

10.REIP紅邊感染點(diǎn)指數(shù)REIP=700+40*[(B4+B7)/2-B5]/(B6-B5)

import numpy as np
# 選擇需要的波段
b4= output.select('B4')# 紅波段
b8= output.select('B8')#近紅外
b2= output.select('B2')#藍(lán)波段
b3= output.select('B3')#綠波段
b7 = output.select('B7')#紅邊3波段
b5= output.select('B5')#紅邊1波段
b6 = output.select('B6')#紅邊2波段# 比值植被指數(shù)RVI = B8/B4
rvi = b8.divide(b4).rename('rvi')
download(rvi)# 紅外植被指數(shù)IPVI = B8/(B8+B4)
ipvi = b8.divide(b8.add(b4)).rename('ipvi')
download(ipvi)# 垂直植被指數(shù)PVI = SIN(45)*B8 - COS(45)*B4
pvi = b8.multiply(np.sin(45)).subtract(b4.multiply(np.cos(45))).rename('pvi')
download(pvi)# 反紅邊葉綠素指數(shù)IRECI=(B8-B4)/(B8+B4)
ireci = b8.subtract(b4).divide(b8.add(b4)).rename('ireci')
download(ireci)# 土壤調(diào)節(jié)植被指數(shù)SAVI=((NIR-R)/(NIR+R+L))(1+L)
L = 0.5  # 土壤校正因子
savi = b8.subtract(b4).divide(b8.add(b4).add(L)).multiply(1 + L).rename('savi')
download(savi)# 大氣阻抗植被指數(shù)ARVI=[NIR - (2xRed-BLUE)] / [NIR + (2xRed-BLUE)]
arvi = b8.subtract(b4.multiply(2).subtract(b2)).divide(b8.add(b4.multiply(2).subtract(b2))).rename('arvi')
download(arvi)# 特定色素簡單比植被指數(shù)PSSRa = B7/B4
pssra = b7.divide(b4).rename('pssra')
download(pssra)# Meris陸地葉綠素指數(shù)MTCI = (B6 - B5)/(B5 - B4 - 0.01)
mtci = b6.subtract(b5).divide(b5.subtract(b4).subtract(0.01)).rename('mtci')
download(mtci)# 修正型葉綠素吸收比值指數(shù)MCARI = (B5 - B4 - 0.2*(B5 - B3))*(B5-B4)
mcari = b5.subtract(b4).subtract((b5.subtract(b3)).multiply(0.2)).multiply(b5.subtract(b4)).rename('mcari')
download(mcari)# REIP紅邊感染點(diǎn)指數(shù)REIP = 700 + 40*((B4+B7)/2 - B5)/(B6 - B5)
reip = ee.Image(700).add(ee.Image(40).multiply(b4.add(b7).divide(2).subtract(b5).divide(b6.subtract(b5)))).rename('reip')
download(reip)

采用gdal計(jì)算各類植被指數(shù)

Sentinel-2 影像文件名 s.tif，然后使用以下命令將指數(shù)計(jì)算轉(zhuǎn)換為 GDAL python本地計(jì)算

安裝gdal庫，安裝gdal庫建議采用本地安裝，去下載whl文件，然后轉(zhuǎn)到放置whl文件的目錄下，pip install 即可
進(jìn)入安裝好gdal庫的虛擬環(huán)境，然后將gdal_calc.py下載地址和s.tif文件放在同一個文件夾下。運(yùn)行python gdal_calc.py文件

RVI：

python gdal_calc.py -A s.tif --A_band=4 -B s_wanzhou.tif --B_band=8 --outfile=rvi_wanzhou.tif --calc="(B/A)"

IPVI：

python gdal_calc.py -A s.tif --A_band=4 -B s_wanzhou.tif --B_band=8 --outfile=ipvi_wanzhou.tif --calc="(B/(A+B))"

PVI：

python gdal_calc.py -A s.tif --A_band=4 -B s.tif --B_band=8 --outfile=pvi.tif --calc="(sin(pi/4)*B)-(cos(pi/4)*A)"

IRECI：

python gdal_calc.py -A s.tif --A_band=4 -B s.tif --B_band=8 --outfile=ireci.tif --calc="((B-A)/(B+A))"

SAVI：

python gdal_calc.py -A s_wanzhou.tif --A_band=4 -B s_wanzhou.tif --B_band=8 --outfile=savi_wanzhou.tif --calc="((B-A)/(B+A+0.5))*(1+0.5)"

ARVI：

python gdal_calc.py -A s_wanzhou.tif -B s_wanzhou.tif -C s_wanzhou.tif --A_band=4 --B_band=8 --C_band=2 --outfile=arvi_wanzhou.tif --calc="((B-(2*A-C))/(B+(2*A-C)))"

PSSRa：

python gdal_calc.py -A s_wanzhou.tif --A_band=4 -B s_wanzhou.tif --B_band=7 --outfile=pssra_wanzhou.tif --calc="(B/A)"

MTCI：

python gdal_calc.py -A s_wanzhou.tif --A_band=4 -B s_wanzhou.tif --B_band=5 -C s_wanzhou.tif --C_band=6 --outfile=mtci_wanzhou.tif --calc="((B-C)/(C-A-0.01))"

MCARI：

python gdal_calc.py -A s_wanzhou.tif --A_band=4 -B s_wanzhou.tif --B_band=5 -C s_wanzhou.tif --C_band=3 --outfile=mcari_wanzhou.tif --calc="((B-A)-(0.2*(B-C)))*(B-A)"

REIP：

python gdal_calc.py -A s_wanzhou.tif --A_band=4 -B s_wanzhou.tif --B_band=5 -C s_wanzhou.tif --C_band=6 -D s_wanzhou.tif --D_band=7 --outfile=reip_wanzhou.tif --calc="(700)+(40*((B+D)/2-A)/(C-A))"

請注意，這里的 s.tif 文件名應(yīng)該與你的實(shí)際文件名相匹配，而且在計(jì)算 SAVI 時，你需要預(yù)先定義 L 的值并將其替換為相應(yīng)的值。

采用gdal_calc.py的命令計(jì)算一下哨兵二號的

NDVI78a (NIR2 – RE3)/ (NIR2 + RE3)

NDVI67 (RE3- RE2)/ (RE3+ RE2)

NDVI58a (NIR2- RE1)/ (NIR2 + RE1)

NDVI56 (RE2- RE1)/ (RE2+ RE1)

NDVI57 (RE3- RE1)/ (RE3+ RE1)

NDVI68a (NIR2 - RE2)/ (NIR2 + RE2)

NDVI48 (NIR - R)/ (NIR + R)
注：藍(lán)波段（B）、綠波段 (G)、紅波段 ?、近紅外波段(NIR1)、窄近紅外波段 (NIR2)、紅邊波段 1 (RE1)、
紅邊波段 2 (RE2)、紅邊波段 3 (RE3).

python gdal_calc.py -A sentinel2.tif --A_band=8 -B sentinel2.tif --B_band=7 --outfile=ndvi78a.tif --calc="(A-B)/(A+B)" --NoDataValue=0python gdal_calc.py -A sentinel2.tif --A_band=7 -B sentinel2.tif --B_band=6 --outfile=ndvi67.tif --calc="(A-B)/(A+B)" --NoDataValue=0python gdal_calc.py -A sentinel2.tif --A_band=8 -B sentinel2.tif --B_band=5 --outfile=ndvi58a.tif --calc="(A-B)/(A+B)" --NoDataValue=0python gdal_calc.py -A sentinel2.tif --A_band=6 -B sentinel2.tif --B_band=5 --outfile=ndvi56.tif --calc="(A-B)/(A+B)" --NoDataValue=0python gdal_calc.py -A sentinel2.tif --A_band=7 -B sentinel2.tif --B_band=5 --outfile=ndvi57.tif --calc="(A-B)/(A+B)" --NoDataValue=0python gdal_calc.py -A sentinel2.tif --A_band=8 -B sentinel2.tif --B_band=6 --outfile=ndvi68a.tif --calc="(A-B)/(A+B)" --NoDataValue=0python gdal_calc.py -A sentinel2.tif --A_band=8 -B sentinel2.tif --B_band=4 --outfile=ndvi48.tif --calc="(A-B)/(A+B)" --NoDataValue=0

1.3 紋理特征參數(shù)提取

采用envi軟件
有研究表明，遙感數(shù)據(jù)的紋理信息能夠增強(qiáng)原始影像亮度的空間信息辨識度，提升地表參數(shù)的反演精度。本研究采用灰度共生矩陣( gray level co-occurrence matrix，GLCM) 提取紋理特征信息，究紋理窗口大小設(shè)置為 3×3，獲取8類 紋 理 特 征 值。
灰度共生矩陣提取紋理特征信息可參考

處理完后可將紋理信息提取出來，可通過APP store 中的Split to Multiple Single-Band Files 工具進(jìn)行直接拆分
Sentinel-2中10個波段每個波段的紋理信息共80個

二.哨兵1號獲取/處理/提取數(shù)據(jù)

哨兵1號數(shù)據(jù)在歐空局中下載，然后采用SNAP軟件對其進(jìn)行一系列處理
可參考鏈接
處理后記得把坐標(biāo)系和投影轉(zhuǎn)成和哨兵2號一樣

2.1 紋理特征參數(shù)提取

同1.3
基于 VH 和VV 極化影像提取紋理特征信息，獲 取 8X2=16 個 紋 理 特 征

三、DEM數(shù)據(jù)

3.1 數(shù)據(jù)下載

可參考這篇文章進(jìn)行數(shù)據(jù)下載

3.2 數(shù)據(jù)處理

一定要注意呀！
獲取的數(shù)據(jù)是30m分辨率的，哨兵數(shù)據(jù)是10米分辨率，在arcGis中搜索resample 需要將DEM重采樣到10米分辨率。
然后在ArcGis中搜索坡度和坡向這兩個工具，分別計(jì)算這兩變量。

四、樣本生物量計(jì)算

代碼如下，d代表胸徑，h代表樹高
查找優(yōu)勢樹種對應(yīng)的異速生長方程寫上就行

import pandas as pd
# 讀取 CSV 文件
df = pd.read_csv('E:\Sentinel12\yangben\yangben.csv', encoding='gbk')
# 定義優(yōu)勢樹種類型及對應(yīng)的
tree_types = {'柏木': lambda d, h: 0.12703 * (d ** 2 * h )** 0.79775,'刺槐': lambda d, h: ( 0.05527 * (d ** 2 * h )** 0.8576) + (  0.02425* (d ** 2 * h )** 0.7908) +( 0.0545 * (d ** 2 * h )** 0.4574),#http://www.xbhp.cn/news/11502.html'櫟類': lambda d, h:0.16625 * ( d ** 2 * h ) ** 0.7821,'柳杉': lambda d, h:( 0.2716 * ( d ** 2 * h ) ** 0.7379 ) + ( 0.0326 * ( d ** 2 * h ) ** 0.8472 ) + ( 0.0250 * ( d ** 2 * h ) ** 1.1778 ) + ( 0.0379 * ( d ** 2 * h ) ** 0.7328 ),'馬尾松': lambda d, h: 0.0973 * (d ** 2 * h )** 0.8285,'其他硬闊': lambda d, h:( 0.0125 * (d ** 2 * h )**1.05 ) + (  0.000933* (d ** 2 * h )**1.23 ) +( 0.000394 * (d ** 2 * h )** 1.20),'杉木': lambda d, h: ( 0.00849 * (d ** 2 * h) ** 1.107230) + ( 0.00175 * (d ** 2 * h )** 1.091916) + 0.00071 * d **3.88664 ,'濕地松': lambda d, h: 0.01016* (d ** 2 * h )**1.098,'香樟': lambda d, h:( 0.05560 * (d ** 2 * h )** 0.850193) + ( 0.00665 * (d ** 2 * h )** 1.051841) +( 0.05987 * (d ** 2 * h )** 0.574327)+( 0.01476 * (d ** 2 * h )** 0.808395) ,'油桐': lambda d, h: ( 0.086217 *d ** 2.00297)+ ( 0.072497 *d ** 2.011502)+( 0.035183 *d ** 1.63929),'楊樹': lambda d, h: ( 0.03 * (d ** 2 * h )** 0.8734) + ( 0.0174 * (d ** 2 * h )** 0.8574) +( 0.4562 * (d ** 2 * h )** 0.3193)+( 0.0028 * (d ** 2 * h )**0.9675 ) ,
}# 遍歷樣本并計(jì)算地上生物量```python
for i, row in df.iterrows():tree_type = row['優(yōu)勢樹']if tree_type in tree_types:sd = row['平均胸徑']h = row['林分平均樹高']biomass = tree_types[tree_type](d, h)df.at[i, '生物量'] = biomass# 將計(jì)算后的地上生物量寫入 CSV 文件
df.to_csv('E:\Sentinel12\yangben\yangben_processed.csv', index=False)

五、樣本變量選取

將樣本導(dǎo)入arcgis,設(shè)置好投影信息，在ArcGis中找到多值提取至點(diǎn)工具。
參考這篇文章

六、隨機(jī)森林建模

參考1
參考2
可以用SPSS進(jìn)行相關(guān)性分析，可參考,選擇相關(guān)性比較大的變量進(jìn)行建模
隨機(jī)森林代碼如下：

6.1導(dǎo)入庫與變量準(zhǔn)備

記得安裝sklearn包
命令行如下：

pip install scikit-learn

本文中設(shè)計(jì)到的所有python代碼均在jupyter notebook中運(yùn)行

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
from pyswarm import pso
import rasterio
# import pydot
import numpy as np
import pandas as pd
import scipy.stats as stats
from pprint import pprint
from sklearn import metrics
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import RandomizedSearchCV
from joblib import dump
# 讀取Excel表格數(shù)據(jù)
data = pd.read_csv(r'E:\Sentinel12\yangben\建模\jianmo.csv' )
y = data.iloc[:, 0].values  # 生物量
X = data.iloc[:, 1:].values  # 指標(biāo)變量
df = pd.DataFrame(data)
random_seed=44
random_forest_seed=np.random.randint(low=1,high=230)
# 分割數(shù)據(jù)集為訓(xùn)練集和測試集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

運(yùn)行完這個代碼塊，可以打印一下，看看變量長什么樣。這里可以看到，選取的變量一共有5個，值分別為1.11691217e+01, 4.37000000e+02, 1.31032691e+00, 7.31299972e-01, 1.13057424e-01 可以打開樣本表看看，這五個變量對應(yīng)的值是否正確。

6.2 選取參數(shù)

決策樹個數(shù)n_estimators

決策樹最大深度max_depth，

最小分離樣本數(shù)（即拆分決策樹節(jié)點(diǎn)所需的最小樣本數(shù)） min_samples_split，

最小葉子節(jié)點(diǎn)樣本數(shù)（即一個葉節(jié)點(diǎn)所需包含的最小樣本數(shù)）min_samples_leaf，

最大分離特征數(shù)（即尋找最佳節(jié)點(diǎn)分割時要考慮的特征變量數(shù)量）max_features

# Search optimal hyperparameter
#對六種超參數(shù)劃定了一個范圍，然后將其分別存入了一個超參數(shù)范圍字典random_forest_hp_range
n_estimators_range=[int(x) for x in np.linspace(start=50,stop=3000,num=60)]
max_features_range=['sqrt','log2',None]
max_depth_range=[int(x) for x in np.linspace(10,500,num=50)]
max_depth_range.append(None)
min_samples_split_range=[2,5,10]
min_samples_leaf_range=[1,2,4,8]random_forest_hp_range={'n_estimators':n_estimators_range,'max_features':max_features_range,'max_depth':max_depth_range,'min_samples_split':min_samples_split_range,'min_samples_leaf':min_samples_leaf_range}
random_forest_model_test_base=RandomForestRegressor()
random_forest_model_test_random=RandomizedSearchCV(estimator=random_forest_model_test_base,param_distributions=random_forest_hp_range,n_iter=200,n_jobs=-1,cv=3,verbose=1,random_state=random_forest_seed)
random_forest_model_test_random.fit(X_train,y_train)best_hp_now=random_forest_model_test_random.best_params_
# Build RF regression model with optimal hyperparameters
random_forest_model_final=random_forest_model_test_random.best_estimator_
print(best_hp_now)

可以看到打印結(jié)果

6.3 誤差分布直方圖

如果直方圖呈現(xiàn)正態(tài)分布或者近似正態(tài)分布，說明模型的預(yù)測誤差分布比較均勻，預(yù)測性能較好。如果直方圖呈現(xiàn)偏斜或者非正態(tài)分布，說明模型在某些情況下的預(yù)測誤差較大，預(yù)測性能可能需要進(jìn)一步改進(jìn)。

# Predict test set data
random_forest_predict=random_forest_model_test_random.predict(X_test)
random_forest_error=random_forest_predict-y_testplt.figure(1)
plt.clf()
plt.hist(random_forest_error)
plt.xlabel('Prediction Error')
plt.ylabel('Count')
plt.grid(False)
plt.show()

6.4 變量重要性可視化展示

# 計(jì)算特征重要性
random_forest_importance = list(random_forest_model_final.feature_importances_)
random_forest_feature_importance = [(feature, round(importance, 8)) for feature, importance in zip(df.columns[4:], random_forest_importance)]
random_forest_feature_importance = sorted(random_forest_feature_importance, key=lambda x:x[1], reverse=True)# 將特征重要性轉(zhuǎn)換為Pyecharts所需的格式
x_data = [item[0] for item in random_forest_feature_importance]
y_data = [item[1] for item in random_forest_feature_importance]# 繪制柱狀圖
bar = (Bar().add_xaxis(x_data).add_yaxis("", y_data).reversal_axis().set_series_opts(label_opts=opts.LabelOpts(position="right")).set_global_opts(title_opts=opts.TitleOpts(title="Variable Importances"),xaxis_opts=opts.AxisOpts(name="Importance",axislabel_opts=opts.LabelOpts(rotate=-45, font_size=10)),yaxis_opts=opts.AxisOpts(name_gap=30, axisline_opts=opts.AxisLineOpts(is_show=False)),tooltip_opts=opts.TooltipOpts(trigger="axis", axis_pointer_type="shadow"))
)bar.render_notebook()

6.5 對精度進(jìn)行驗(yàn)證

這里可輸出相關(guān)的精度值

# Verify the accuracy
random_forest_pearson_r=stats.pearsonr(y_test,random_forest_predict)
random_forest_R2=metrics.r2_score(y_test,random_forest_predict)
random_forest_RMSE=metrics.mean_squared_error(y_test,random_forest_predict)**0.5
print('Pearson correlation coefficient = {0}、R2 = {1} and RMSE = {2}.'.format(random_forest_pearson_r[0],random_forest_R2,random_forest_RMSE))

6.6 預(yù)測生物量

注意幾個tif數(shù)據(jù)的nodata值不一樣，最好全部都將nodata值設(shè)為0，最后得到的biomass圖像按照掩膜提取，nodata就變回來啦

import numpy as np
import rasterio
from tqdm import tqdm# 讀取五個柵格指標(biāo)變量數(shù)據(jù)并整合為一個矩陣
with rasterio.open(r'E:\Sentinel12\yangben\建模\result\nodata\podu.tif') as src:data1 = src.read(1)meta = src.meta
with rasterio.open(r'E:\Sentinel12\yangben\建模\result\nodata\dem.tif') as src:data2 = src.read(1)
with rasterio.open(r'E:\Sentinel12\yangben\建模\result\nodata\lai.tif') as src:data3 = src.read(1)
with rasterio.open(r'E:\Sentinel12\yangben\建模\result\nodata\ndvi.tif') as src:data4 = src.read(1)
with rasterio.open(r'E:\Sentinel12\yangben\建模\result\nodata\pvi.tif') as src:data5 = src.read(1)X = np.stack((data1, data2, data3, data4, data5), axis=-1)# 清洗輸入數(shù)據(jù)
X_2d = X.reshape(-1, X.shape[-1])
# 檢查數(shù)據(jù)中是否存在 NaN 值
print(np.isnan(X_2d).any())  # True# 將 nodata 值替換為0
X_2d[np.isnan(X_2d)] = 0# 使用已經(jīng)訓(xùn)練好的隨機(jī)森林模型對整合后的柵格指標(biāo)變量數(shù)據(jù)進(jìn)行生物量的預(yù)測
y_pred = []
for i in tqdm(range(0, X_2d.shape[0], 10000)):y_pred_chunk = random_forest_model_test_random.predict(X_2d[i:i+10000])y_pred.append(y_pred_chunk)
y_pred = np.concatenate(y_pred)# 將預(yù)測結(jié)果保存為一個新的柵格數(shù)據(jù)
with rasterio.open(r'E:\Sentinel12\yangben\建模\result\biomass_v5.tif', 'w', **meta) as dst:dst.write(y_pred.reshape(X.shape[:-1]), 1)
print("預(yù)測結(jié)束")

運(yùn)行之后會在下面出現(xiàn)進(jìn)度條，進(jìn)度條走完了就代碼預(yù)測完了

七、生物量制圖

7.1. 將得到的biomass.tif文件按掩膜提取

 掩膜文件選擇用于預(yù)測的tif文件，目的是將nodata值還原回來。

7.2. 提取森林掩膜區(qū)域

[參考鏈接](https://www.bilibili.com/video/BV1qv4y1A75B?p=10&vd_source=ddb0eb9d8fde0d0629fc9f25dc6915e5) 

• 這里需要全球土地覆蓋10米的更精細(xì)分辨率觀測和監(jiān)測（FROM-GLC10）數(shù)據(jù)。

• 我們在GEE平臺上下載研究區(qū)的GLC10圖，參考鏈接

1. 在 Google Earth Engine Code Editor 中搜索“ESA Global Land Cover 10m v2.0.0”并加載該數(shù)據(jù)集。

var glc = ee.Image('ESA/GLOBCOVER_L4_200901_200912_V2_3');

2. 定義您感興趣的區(qū)域。
   這里的table需要先將shp文件上傳到平臺，再點(diǎn)那個小箭頭導(dǎo)入
   
   這里就會有table出現(xiàn)
   

var geometry = table;

3. 使用 Image.clip() 函數(shù)將圖像裁剪為您感興趣的區(qū)域。

var clipped = glc.clip(geometry);

4. 使用 Export.image.toDrive() 函數(shù)導(dǎo)出圖像。以下代碼將圖像導(dǎo)出到您的 Google Drive 中。

Export.image.toDrive({image: clipped,description: 'GLC',folder: 'my_folder',scale: 10,region: geometry
});

其中，description 是導(dǎo)出圖像的名稱，folder 是導(dǎo)出圖像的文件夾，scale 是導(dǎo)出圖像的分辨率，region 是導(dǎo)出圖像的區(qū)域。

5. 點(diǎn)擊“Run”按鈕運(yùn)行代碼。代碼運(yùn)行完成后，您可以在 Google Drive 中找到導(dǎo)出的圖像文件。

• 20序號代表森林，按屬性提取，可以把森林提取出來，按屬性提取工具。
  
  將生物量.tif按照掩膜tif掩膜即可得到森林區(qū)域的生物量。

八. 需要注意的點(diǎn)

1. 每個柵格變量數(shù)據(jù)一定要保持行數(shù)和列數(shù)一致，這不僅是為了”多值提取至點(diǎn)“等一一對應(yīng)，更是為了最后估算生物量的時候生成二維矩陣輸入模型，保證輸入數(shù)據(jù)的維度一致。
2. 第一步的前提是第二部，投影！投影！投影!重要的事情說三遍，一定要在相同坐標(biāo)系下面
3. 注意nodata值的處理，最好所有的影像nodata設(shè)置為0,這樣最后輸入模型中可以減少計(jì)算量。