第一次提交
This commit is contained in:
tangchao
164
.gitignore
vendored
Normal file
164
.gitignore
vendored
Normal file
@ -0,0 +1,164 @@
|
||||
# 唐超添加
|
||||
/.idea
|
||||
|
||||
|
||||
# Byte-compiled / optimized / DLL files
|
||||
__pycache__/
|
||||
*.py[cod]
|
||||
*$py.class
|
||||
|
||||
# C extensions
|
||||
*.so
|
||||
|
||||
# Distribution / packaging
|
||||
.Python
|
||||
build/
|
||||
develop-eggs/
|
||||
dist/
|
||||
downloads/
|
||||
eggs/
|
||||
.eggs/
|
||||
lib/
|
||||
lib64/
|
||||
parts/
|
||||
sdist/
|
||||
var/
|
||||
wheels/
|
||||
share/python-wheels/
|
||||
*.egg-info/
|
||||
.installed.cfg
|
||||
*.egg
|
||||
MANIFEST
|
||||
|
||||
# PyInstaller
|
||||
# Usually these files are written by a python script from a template
|
||||
# before PyInstaller builds the exe, so as to inject date/other infos into it.
|
||||
*.manifest
|
||||
*.spec
|
||||
|
||||
# Installer logs
|
||||
pip-log.txt
|
||||
pip-delete-this-directory.txt
|
||||
|
||||
# Unit test / coverage reports
|
||||
htmlcov/
|
||||
.tox/
|
||||
.nox/
|
||||
.coverage
|
||||
.coverage.*
|
||||
.cache
|
||||
nosetests.xml
|
||||
coverage.xml
|
||||
*.cover
|
||||
*.py,cover
|
||||
.hypothesis/
|
||||
.pytest_cache/
|
||||
cover/
|
||||
|
||||
# Translations
|
||||
*.mo
|
||||
*.pot
|
||||
|
||||
# Django stuff:
|
||||
*.log
|
||||
local_settings.py
|
||||
db.sqlite3
|
||||
db.sqlite3-journal
|
||||
|
||||
# Flask stuff:
|
||||
instance/
|
||||
.webassets-cache
|
||||
|
||||
# Scrapy stuff:
|
||||
.scrapy
|
||||
|
||||
# Sphinx documentation
|
||||
docs/_build/
|
||||
|
||||
# PyBuilder
|
||||
.pybuilder/
|
||||
target/
|
||||
|
||||
# Jupyter Notebook
|
||||
.ipynb_checkpoints
|
||||
|
||||
# IPython
|
||||
profile_default/
|
||||
ipython_config.py
|
||||
|
||||
# pyenv
|
||||
# For a library or package, you might want to ignore these files since the code is
|
||||
# intended to run in multiple environments; otherwise, check them in:
|
||||
# .python-version
|
||||
|
||||
# pipenv
|
||||
# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
|
||||
# However, in case of collaboration, if having platform-specific dependencies or dependencies
|
||||
# having no cross-platform support, pipenv may install dependencies that don't work, or not
|
||||
# install all needed dependencies.
|
||||
#Pipfile.lock
|
||||
|
||||
# poetry
|
||||
# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
|
||||
# This is especially recommended for binary packages to ensure reproducibility, and is more
|
||||
# commonly ignored for libraries.
|
||||
# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
|
||||
#poetry.lock
|
||||
|
||||
# pdm
|
||||
# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
|
||||
#pdm.lock
|
||||
# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
|
||||
# in version control.
|
||||
# https://pdm.fming.dev/#use-with-ide
|
||||
.pdm.toml
|
||||
|
||||
# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
|
||||
__pypackages__/
|
||||
|
||||
# Celery stuff
|
||||
celerybeat-schedule
|
||||
celerybeat.pid
|
||||
|
||||
# SageMath parsed files
|
||||
*.sage.py
|
||||
|
||||
# Environments
|
||||
.env
|
||||
.venv
|
||||
env/
|
||||
venv/
|
||||
ENV/
|
||||
env.bak/
|
||||
venv.bak/
|
||||
|
||||
# Spyder project settings
|
||||
.spyderproject
|
||||
.spyproject
|
||||
|
||||
# Rope project settings
|
||||
.ropeproject
|
||||
|
||||
# mkdocs documentation
|
||||
/site
|
||||
|
||||
# mypy
|
||||
.mypy_cache/
|
||||
.dmypy.json
|
||||
dmypy.json
|
||||
|
||||
# Pyre type checker
|
||||
.pyre/
|
||||
|
||||
# pytype static type analyzer
|
||||
.pytype/
|
||||
|
||||
# Cython debug symbols
|
||||
cython_debug/
|
||||
|
||||
# PyCharm
|
||||
# JetBrains specific template is maintained in a separate JetBrains.gitignore that can
|
||||
# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
|
||||
# and can be added to the global gitignore or merged into this file. For a more nuclear
|
||||
# option (not recommended) you can uncomment the following to ignore the entire idea folder.
|
||||
#.idea/
|
||||
72
Setup_Window.m
Normal file
72
Setup_Window.m
Normal file
@ -0,0 +1,72 @@
|
||||
function [windows] = Setup_Window()
|
||||
|
||||
windows = [];
|
||||
% -----Attention: fitting window range for sFLD<EFBFBD><EFBFBD>3FLD<EFBFBD><EFBFBD>SFM<EFBFBD><EFBFBD>SVD<EFBFBD><EFBFBD>DOAS;----- %
|
||||
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>sFLD, iFLD, pFLD
|
||||
FLD = [];
|
||||
FLD.outerH = [656.0, 656.2]; % <EFBFBD>˵㲨<EFBFBD>Σ<EFBFBD>ȡ<EFBFBD><EFBFBD>ֵ
|
||||
FLD.outerB = [687.0, 687.2];
|
||||
FLD.outerW = [717.0, 717.2];
|
||||
FLD.outerA = [758.6, 758.8];
|
||||
FLD.windowH = [656, 659]; % <EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сֵ
|
||||
FLD.windowB = [686, 689];
|
||||
FLD.windowW = [717, 722];
|
||||
FLD.windowA = [758, 770];
|
||||
windows.FLD = FLD;
|
||||
|
||||
% 3FLD
|
||||
FLD3 = [];
|
||||
FLD3.FLD3_leftH = [656.0, 656.2]; FLD3.FLD3_rightH = [659.0, 659.2];
|
||||
FLD3.FLD3_leftB = [687.0, 687.2]; FLD3.FLD3_rightB = [689.0, 689.2];
|
||||
FLD3.FLD3_leftW = [717.0, 717.2]; FLD3.FLD3_rightW = [722.0, 722.2];
|
||||
FLD3.FLD3_leftA = [758.6, 758.8]; FLD3.FLD3_rightA = [770.0, 770.2];
|
||||
FLD3.windowH = [656, 659]; % <EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сֵ
|
||||
FLD3.windowB = [686, 689];
|
||||
FLD3.windowW = [717, 722];
|
||||
FLD3.windowA = [758, 770];
|
||||
windows.FLD3 = FLD3;
|
||||
|
||||
% SFM
|
||||
SFM = [];
|
||||
SFM.SFM_windowH = [656, 667]; % <EFBFBD><EFBFBD><EFBFBD>η<EFBFBD>Χ
|
||||
SFM.SFM_windowB = [687, 698];
|
||||
SFM.SFM_windowW = [717, 728];
|
||||
SFM.SFM_windowA = [759, 770];
|
||||
SFM.SFM_Ha = [656.0 656.2]; % ij<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ӫ<EFBFBD><EFBFBD>ֵ,SFM<EFBFBD><EFBFBD>SVD<EFBFBD>㷨ֱ<EFBFBD>Ӷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||||
SFM.SFM_O2B = [687.0 687.2];
|
||||
SFM.SFM_H2O = [720.0 720.2];
|
||||
SFM.SFM_O2A = [760.0 760.2];
|
||||
windows.SFM = SFM;
|
||||
% SVD
|
||||
SVD = [];
|
||||
SVD.SVD_windowH = [655, 670];
|
||||
SVD.SVD_windowB = [675, 710]; % wide O2B and O2A fitting window; % e.g. narrow O2A fitting window:[759.86 762.79];
|
||||
SVD.SVD_windowW = [700, 740];
|
||||
SVD.SVD_windowA = [740, 780];
|
||||
SVD.SVD_Ha = [656.0 656.2]; % ij<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ӫ<EFBFBD><EFBFBD>ֵ,SFM<EFBFBD><EFBFBD>SVD<EFBFBD>㷨ֱ<EFBFBD>Ӷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||||
SVD.SVD_O2B = [687.0 687.2];
|
||||
SVD.SVD_H2O = [720.0 720.2];
|
||||
SVD.SVD_O2A = [760.0 760.2];
|
||||
% Example SVD settings
|
||||
% if rolling > 0, abs(rolling) = the number of spectra in a centered moving window used for training the SVD
|
||||
% if rolling = 0, it will use all spectra from the day for training the SVD
|
||||
% if rolling < 0, it will use that number of spectra but evenly dispersed throughout the day for training the SVD
|
||||
SVD.rolling = 0; % moving window of 5 spectra
|
||||
windows.SVD = SVD;
|
||||
|
||||
% DOAS, F_SFM, <EFBFBD><EFBFBD><EFBFBD>ɷַ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD>η<EFBFBD>Χ<EFBFBD>ڷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сֵ
|
||||
Others = []; % <EFBFBD><EFBFBD><EFBFBD>ַ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>õIJ<EFBFBD><EFBFBD>η<EFBFBD>Χ<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD>Сֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||||
Others.WindowH = [656, 659]; % Attention: ע<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||||
Others.WindowB = [686, 689];
|
||||
Others.WindowW = [717, 722];
|
||||
Others.WindowA = [758, 770];
|
||||
windows.Others = Others;
|
||||
|
||||
% -----Attention: Modify the range of spectrum;----- %
|
||||
% --<EFBFBD>õ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ⱥͷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
|
||||
windows.wl_680 = [676, 680]; % 676-680nm, red
|
||||
windows.wl_705 = [703, 707]; % 703-707nm, red
|
||||
windows.wl_740 = [736, 740]; % 736-740nm, rededge
|
||||
windows.wl_770 = [770, 774]; % 770-774nm, near_red
|
||||
windows.wl_800 = [795, 799]; % 795-799nm, near_red
|
||||
end
|
||||
310
sif_methods.py
Normal file
310
sif_methods.py
Normal file
@ -0,0 +1,310 @@
|
||||
import numpy as np
|
||||
from scipy.special import eval_legendre
|
||||
from scipy import optimize
|
||||
# 所有提取方法均基于xarray数据，xarray数据的变量与波长、时间、光谱仪绑定
|
||||
def cal_inside_bands_ave(data):
|
||||
'''
|
||||
根据多个光谱找出窗口内数据最低点对应波长
|
||||
'''
|
||||
sky_spec = data.sky
|
||||
veg_spec = data.veg
|
||||
wvl_inside_band_l = np.mean(veg_spec.idxmin(dim='Wavelength')).values
|
||||
wvl_inside_band_e = np.mean(sky_spec.idxmin(dim='Wavelength')).values
|
||||
return[wvl_inside_band_l,wvl_inside_band_e]
|
||||
|
||||
def cal_outside_values_mean(data,outer):
|
||||
'''
|
||||
计算肩部窗口的均值
|
||||
'''
|
||||
_data = data.where((data.Wavelength>outer[0])&(data.Wavelength<outer[1]),drop=True)
|
||||
wvl_outer_mean = _data['Wavelength'].mean().values
|
||||
_mean = _data.mean(dim = 'Wavelength')
|
||||
E_out = _mean.sky.values
|
||||
L_out = _mean.veg.values
|
||||
return L_out,E_out,wvl_outer_mean
|
||||
|
||||
def sfld(data,wl_range,outer):
|
||||
|
||||
"""
|
||||
Standard FLD (Huaize Feng)
|
||||
|
||||
input:
|
||||
data,xarray dataset
|
||||
wl_range, a window around Fraunhofer lines position, [start,end]
|
||||
outer, a window outside the Fraunhofer lines, [start,end]
|
||||
"""
|
||||
sif = []
|
||||
refl = []
|
||||
nmeas_ = data.Measures.size
|
||||
|
||||
data = data.where((data.Wavelength>wl_range[0])&(data.Wavelength<wl_range[1]),drop=True)
|
||||
[wvl_inside_band_l,wvl_inside_band_e]=cal_inside_bands_ave(data)
|
||||
for i in range(0,nmeas_):
|
||||
_data = data.isel(Measures=i)
|
||||
veg_out,sky_out,_ = cal_outside_values_mean(_data,outer)
|
||||
veg_in = _data.veg.sel(Wavelength = wvl_inside_band_l,method='nearest').values
|
||||
sky_in = _data.sky.sel(Wavelength = wvl_inside_band_e,method='nearest').values
|
||||
_sif = (sky_out*veg_in - sky_in*veg_out)/(sky_out - sky_in)
|
||||
_refl = (veg_in - _sif)/ sky_in
|
||||
sif.append(_sif)
|
||||
refl.append(_refl)
|
||||
return[sif,refl]
|
||||
|
||||
def fld3(data,wl_range,outer_left,outer_right):
|
||||
"""
|
||||
3FLD (Huaize Feng)
|
||||
input:
|
||||
data,xarray dataset
|
||||
wl_range, a window around Fraunhofer lines position, [start,end]
|
||||
outer_*, a window outside the Fraunhofer lines, [start,end]
|
||||
outer_left, left window outside the Fraunhofer lines
|
||||
outer_right, left window outside the Fraunhofer lines
|
||||
return:
|
||||
list[sif,reflectance]
|
||||
"""
|
||||
sif = []
|
||||
refl = []
|
||||
nmeas_ = data.Measures.size
|
||||
|
||||
data = data.where((data.Wavelength>wl_range[0])&(data.Wavelength<wl_range[1]),drop=True)
|
||||
[wvl_inside_band_l,wvl_inside_band_e]=cal_inside_bands_ave(data)
|
||||
for i in range(0,nmeas_):
|
||||
_data = data.isel(Measures=i)
|
||||
veg_out_left,sky_out_left,wvl_outer_left = cal_outside_values_mean(_data,outer_left)
|
||||
veg_out_right,sky_out_right,wvl_outer_right = cal_outside_values_mean(_data,outer_right)
|
||||
veg_in = _data.veg.sel(Wavelength = wvl_inside_band_l,method='nearest').values
|
||||
sky_in = _data.sky.sel(Wavelength = wvl_inside_band_e,method='nearest').values
|
||||
|
||||
# 根据离吸收峰的距离的反比进行赋权
|
||||
wight_left = (wvl_outer_right - wvl_inside_band_e)/(wvl_outer_right - wvl_outer_left)
|
||||
wight_right = (wvl_inside_band_e - wvl_outer_left)/(wvl_outer_right - wvl_outer_left)
|
||||
|
||||
_sif = (veg_in - (sky_in/((wight_left*sky_out_left) + (wight_right*sky_out_right))) * ((wight_left*veg_out_left) + (wight_right*veg_out_right))) / (1-(sky_in / ((wight_left*sky_out_left) + (wight_right*sky_out_right))))
|
||||
_refl = (veg_in - _sif)/ sky_in
|
||||
sif.append(_sif)
|
||||
refl.append(_refl)
|
||||
return[sif,refl]
|
||||
|
||||
def sfm(data,wl_range,band):
|
||||
|
||||
"""
|
||||
Spectral Fitting Method (Huaize Feng)
|
||||
input:
|
||||
data,
|
||||
xarray dataset
|
||||
wl_range,[start,end]
|
||||
a window around Fraunhofer lines position, ,about 10 nm
|
||||
for fitting a ployniam or gussian function.
|
||||
band,float/int
|
||||
exact position of the Fraunhofer line
|
||||
return:
|
||||
list[sif,reflectance,rmse,B]
|
||||
sif, numpy array, sif at the Fraunhofer line for each measurement
|
||||
reflectance, numpy array
|
||||
rmse, RMSE
|
||||
B, numpy array, the parameters of the fitting equation
|
||||
[a,b,c,d,e,f]: a, b, c for polynimal refelectance
|
||||
d, float, MAX sif value [0,10] mw/m2/nm/sr
|
||||
e, position of MAX sif value - wl_range, [0,wavelength.size], nm
|
||||
f, full width at half maximum, [0,wavelength.size*5]
|
||||
"""
|
||||
sif,refl,rmse,B = [],[],[],[]
|
||||
data = data.where((data.Wavelength>wl_range[0])&(data.Wavelength<wl_range[1]),drop=True)
|
||||
abosorb_line_position = np.where(data.Wavelength == data.Wavelength.sel(Wavelength=760,method='nearest').values)
|
||||
_nmeas = data.Measures.size
|
||||
_nwvl = data.Wavelength.size
|
||||
_x = (data.Wavelength -np.min(data.Wavelength)).values
|
||||
poly_refl = [_x**2,_x,np.ones(_nwvl)]
|
||||
poly_sif = [_x**2,_x,np.ones(_nwvl)]
|
||||
|
||||
for i in range(0,_nmeas):
|
||||
sky_spec = data.sky[i].values
|
||||
veg_spec = data.veg[i].values
|
||||
_X = np.concatenate([poly_refl*sky_spec,poly_refl]).T
|
||||
_B = np.linalg.lstsq(_X,veg_spec,rcond=-1)
|
||||
_sif = np.array(poly_sif).T.dot(_B[0][-3:])
|
||||
_refl = (veg_spec - _sif) / (sky_spec+0.000001)
|
||||
sif.append(_sif[abosorb_line_position][0])
|
||||
refl.append(_refl[abosorb_line_position][0])
|
||||
rmse.append(np.sqrt(np.sum(_B[1]**2)/_nwvl))
|
||||
B.append(_B[0])
|
||||
return [sif,refl,rmse,B]
|
||||
|
||||
|
||||
def f(x, a, b, c, d, e, f):
|
||||
sky_spec = x['sky_spec']
|
||||
_x = x['_x']
|
||||
refl = a * _x**2*sky_spec + b* _x*sky_spec + c*sky_spec
|
||||
sif = d*np.exp(-(_x-e)**2/f)
|
||||
y_hat = refl + sif
|
||||
return y_hat
|
||||
|
||||
def f_cal(X, a,b,c,d,e,f):
|
||||
sky_spec = X['sky_spec']
|
||||
_x = X['_x']
|
||||
refl = a * _x**2*sky_spec + b* _x*sky_spec + c*sky_spec
|
||||
sif = d*np.exp(-(_x-e)**2/f)
|
||||
y_hat = refl + sif
|
||||
return sif,refl
|
||||
|
||||
def sfm_gaussian(data,wl_range,band=760):
|
||||
|
||||
'''
|
||||
Spectral Fitting Method (Gaussian Ver. (Huaize Feng))
|
||||
input:
|
||||
data,
|
||||
xarray dataset
|
||||
wl_range,[start,end]
|
||||
a window around Fraunhofer lines position, ,about 10 nm
|
||||
for fitting a polynomial or gaussian function.
|
||||
band,float/int
|
||||
exact position of the Fraunhofer line
|
||||
return:
|
||||
list[sif,reflectance,rmse,B]
|
||||
sif, numpy array, sif at the Fraunhofer line for each measurement
|
||||
reflectance, numpy array
|
||||
rmse, RMSE
|
||||
B, numpy array, the parameters of the fitting equation
|
||||
[a,b,c,d,e,f]: a, b, c for polynimal refelectance
|
||||
d, float, MAX sif value [0,10] mw/m2/nm/sr
|
||||
e, position of MAX sif value - wl_range, [0,wavelength.size], nm
|
||||
f, full width at half maximum
|
||||
'''
|
||||
sif,refl,rmse,B = [],[],[],[]
|
||||
data = data.where((data.Wavelength>wl_range[0])&(data.Wavelength<wl_range[1]),drop=True)
|
||||
abosorb_line_position = np.where(data.Wavelength == data.Wavelength.sel(Wavelength=band,method='nearest').values)
|
||||
_nmeas = data.Measures.size
|
||||
_nwvl = data.Wavelength.size
|
||||
_x = (data.Wavelength -np.min(data.Wavelength)).values
|
||||
poly_refl = [_x**2,_x,np.ones(_nwvl)]
|
||||
poly_sif = [_x**2,_x,np.ones(_nwvl)]
|
||||
|
||||
bounds=((-np.inf, -np.inf,-np.inf, 0, 0, 0), (np.inf, np.inf,np.inf, 10, _nwvl, _nwvl*5))
|
||||
|
||||
for i in range(0,_nmeas):
|
||||
sky_spec = data.sky[i].values
|
||||
veg_spec = data.veg[i].values
|
||||
_X = {'sky_spec':sky_spec,'_x':_x}
|
||||
_B = optimize.curve_fit(f, _X, veg_spec,bounds=bounds)[0]
|
||||
# print(_B)
|
||||
|
||||
_sif,_ = f_cal(_X,*_B)
|
||||
# print(_sif)
|
||||
_rmse = np.sqrt(np.sum((f(_X,*_B)-veg_spec)**2)/_nwvl)
|
||||
_refl = (veg_spec - _sif) / sky_spec
|
||||
sif.append(_sif[abosorb_line_position][0])
|
||||
refl.append(_refl[abosorb_line_position][0])
|
||||
rmse.append( _rmse)
|
||||
B.append(_B[0])
|
||||
return [sif,refl,rmse,B]
|
||||
|
||||
def doas(data,wl_range,band=760):
|
||||
|
||||
"""
|
||||
Spectral Fitting Method (Huaize Feng)
|
||||
input:
|
||||
data,
|
||||
xarray dataset
|
||||
wl_range,[start,end]
|
||||
a window around Fraunhofer lines position, ,about 10 nm
|
||||
for fitting a ployniam or gussian function.
|
||||
band,float/int
|
||||
exact position of the Fraunhofer line
|
||||
return:
|
||||
list[sif,reflectance,rmse,B]
|
||||
sif, numpy array, sif at the Fraunhofer line for each measurement
|
||||
reflectance, numpy array
|
||||
rmse, RMSE
|
||||
B, numpy array, the parameters of the fitting equation
|
||||
[a,b,c,d,e,f]: a, b, c for polynimal refelectance
|
||||
d, float, MAX sif value [0,10] mw/m2/nm/sr
|
||||
e, position of MAX sif value - wl_range, [0,wavelength.size], nm
|
||||
f, full width at half maximum, [0,wavelength.size*5]
|
||||
"""
|
||||
sif,refl,rmse,B = [],[],[],[]
|
||||
data = data.where((data.Wavelength>wl_range[0])&(data.Wavelength<wl_range[1]),drop=True)
|
||||
absorb_line_position = np.where(data.Wavelength == data.Wavelength.sel(Wavelength=band,method='nearest').values)
|
||||
_nmeas = data.Measures.size
|
||||
_wvl = data.Wavelength.values
|
||||
_nwvl = _wvl.size
|
||||
_x = np.interp(_wvl, (_wvl.min(),_wvl.max()),(-1,1))
|
||||
# normalize standard SIF template by the mean value
|
||||
_hf = data.hf.values/(np.mean(data.hf.values))
|
||||
# create base function by legendre polynomial equations
|
||||
_ref_base = np.array([eval_legendre(n, _x) for n in np.arange(6)])
|
||||
|
||||
for i in range(0,_nmeas):
|
||||
print(data.Measures[i].values)
|
||||
sky_spec = data.sky[i].values
|
||||
veg_spec = data.veg[i].values
|
||||
_X = np.concatenate([_ref_base,(_hf/(veg_spec+0.000001111)).reshape(1,-1)]).T
|
||||
_y = np.log(veg_spec+0.000001111) - np.log(sky_spec+0.00000111111)
|
||||
_B = np.linalg.lstsq(_X,_y,rcond=-1)
|
||||
_sif = _hf * _B[0][-1]
|
||||
_refl = (veg_spec - _sif) / (sky_spec+0.0000011111)
|
||||
sif.append(_sif[absorb_line_position][0])
|
||||
refl.append(_refl[absorb_line_position][0])
|
||||
rmse.append(np.sqrt(np.sum(_B[1]**2)/_nwvl))
|
||||
B.append(_B[0])
|
||||
|
||||
return [sif,refl,rmse,B]
|
||||
|
||||
def svd(data,wl_range,band=760,num_vector=20,pow_of_refl = 5):
|
||||
|
||||
"""
|
||||
Singular Vector Decomposition (Huaize Feng)
|
||||
input:
|
||||
data,
|
||||
xarray dataset
|
||||
wl_range,[start,end]
|
||||
a window around Fraunhofer lines position, ,about 10 nm
|
||||
for fitting a polynomial or gaussian function.
|
||||
band,float/int
|
||||
exact position of the Fraunhofer line
|
||||
num_vector,int
|
||||
pca保留的观测维度（个数），相当于去除大气的噪音
|
||||
pow_of_refl,int
|
||||
反射率所用多项式的最高次数
|
||||
Highest power of the polynomial equation denoting reflectance
|
||||
return:
|
||||
list[sif,reflectance,rmse,B]
|
||||
sif, numpy array, sif at the Fraunhofer line for each measurement
|
||||
reflectance, numpy array
|
||||
rmse, RMSE
|
||||
B, numpy array, the parameters of the fitting equation
|
||||
[a,b,c,d,e,f]: a, b, c for polynimal refelectance
|
||||
d, float, MAX sif value [0,10] mw/m2/nm/sr
|
||||
e, position of MAX sif value - wl_range, [0,wavelength.size], nm
|
||||
f, full width at half maximum, [0,wavelength.size*5]
|
||||
"""
|
||||
sif,refl,rmse,B = [],[],[],[]
|
||||
data = data.where((data.Wavelength>wl_range[0])&(data.Wavelength<wl_range[1]),drop=True)
|
||||
abosorb_line_position = np.where(data.Wavelength == data.Wavelength.sel(Wavelength=band,method='nearest').values)
|
||||
_nmeas = data.Measures.size
|
||||
_wvl = data.Wavelength.values
|
||||
_nwvl = _wvl.size
|
||||
_hf = data.hf.values # 重采样后的标准sif
|
||||
|
||||
u, s, vh = np.linalg.svd(data.sky) # Svd分解----------------------------------------------------------------------------------------------------
|
||||
v = -vh
|
||||
|
||||
tmp1=(_wvl-np.mean(_wvl)).reshape(-1,1)
|
||||
tmp2=np.arange(pow_of_refl+1)
|
||||
p1 = np.power(tmp1, tmp2)*v[0].reshape(-1,1)
|
||||
p2 = np.power((_wvl-np.mean(_wvl)).reshape(-1,1), np.arange(pow_of_refl+1))*v[1].reshape(-1,1)
|
||||
p3 = v[:,2:num_vector]
|
||||
_X = np.concatenate([p1,p2,p3,_hf.reshape(-1,1)],axis=1)
|
||||
|
||||
for i in range(0,_nmeas):
|
||||
sky_spec = data.sky[i].values
|
||||
veg_spec = data.veg[i].values
|
||||
_y = veg_spec
|
||||
_B = np.linalg.lstsq(_X,_y,rcond=-1) # https://www.zhihu.com/question/40540185?sort=created
|
||||
_sif = _hf * _B[0][-1]
|
||||
_refl = (veg_spec - _sif) / (sky_spec + 0.0000001)
|
||||
sif.append(_sif[abosorb_line_position][0])
|
||||
refl.append(_refl[abosorb_line_position][0])
|
||||
rmse.append(np.sqrt(np.sum(_B[1]**2)/_nwvl))
|
||||
B.append(_B[0])
|
||||
return [sif,refl,rmse,B]
|
||||
|
||||
251
sif_retrieval.py
Normal file
251
sif_retrieval.py
Normal file
@ -0,0 +1,251 @@
|
||||
# Multi Processing Version.
|
||||
# 2022/1/14 14:52
|
||||
|
||||
# 唐超：如何将此文件打包成一个exe
|
||||
# （1）cd D:\csharp_vs2017\easySif\sifAlgorithm_Python_exe
|
||||
# （2）D:\software\Anaconda3\envs\python_36\Scripts\pyi-makespec.exe -F D:\PycharmProjects\sif\sif_retrieval.py
|
||||
# （3）在生成的spec文件中加入import sys
|
||||
# sys.setrecursionlimit(5000)
|
||||
# （4）D:\software\Anaconda3\envs\python_36\Scripts\pyinstaller.exe D:\csharp_vs2017\easySif\sifAlgorithm_Python_exe\sif_retrieval.spec
|
||||
|
||||
import os
|
||||
import sys
|
||||
from glob import glob
|
||||
import ast
|
||||
import csv
|
||||
import sif_methods
|
||||
import argparse
|
||||
from multiprocessing import Pool
|
||||
import multiprocessing
|
||||
from scipy.interpolate import splev, splrep
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
try:
|
||||
import xarray as xr
|
||||
except Exception as e:
|
||||
print(e)
|
||||
cmd = 'pip install xarray -i https://pypi.tuna.tsinghua.edu.cn/simple --user'
|
||||
os.system(cmd)
|
||||
import xarray as xr
|
||||
finally:
|
||||
pass
|
||||
|
||||
|
||||
def intep(standard, x2):
|
||||
'''
|
||||
插值，根据给定的波长序列将标准SIF重新插值
|
||||
输入：标准SIF，两列，一列波长一列值；给定新的波长序列
|
||||
输出：插值后的SIF值
|
||||
'''
|
||||
spl = splrep(standard.index, standard)
|
||||
y2 = splev(x2, spl)
|
||||
return y2
|
||||
|
||||
|
||||
def sel_method(x):
|
||||
'''
|
||||
选择方法，
|
||||
输入：方法字符串
|
||||
输出：方法
|
||||
'''
|
||||
return {
|
||||
'sfld': sif_methods.sfld,
|
||||
'3fld': sif_methods.fld3,
|
||||
'svd': sif_methods.svd,
|
||||
'sfm': sif_methods.sfm,
|
||||
'sfm_gaussian': sif_methods.sfm_gaussian,
|
||||
'doas': sif_methods.doas
|
||||
}.get(x, sif_methods.sfld)
|
||||
|
||||
|
||||
def read_files(f):
|
||||
'''
|
||||
使用csv.reader读取文件
|
||||
'''
|
||||
with open(f) as csvfile:
|
||||
data = list(csv.reader(csvfile)) # 读csv并转化为list
|
||||
# data = np.array(data)
|
||||
data = np.array(data, dtype=object) # 转化为np数组
|
||||
data = [np.array(d) for d in data]
|
||||
csvfile.close()
|
||||
return data
|
||||
|
||||
|
||||
def Convert(lst):
|
||||
'''
|
||||
针对<键名 键值...>组织形式的文件头解析
|
||||
'''
|
||||
# tc
|
||||
# res_dct = {lst[i]: lst[i + 1] for i in range(0, lst.size-1, 2)}
|
||||
res_dct = {lst[i]: lst[i + 1] for i in range(0, len(lst)-1, 2)}
|
||||
return res_dct
|
||||
|
||||
# parsing
|
||||
|
||||
|
||||
def data_parser(data, dt):
|
||||
shp = np.arange(len(data)).tolist()
|
||||
_iter = iter(shp)
|
||||
|
||||
# 读取csv的元信息（第一行）
|
||||
data[next(_iter)][1:] # 略过csv的第一行
|
||||
header_dict = {'Measure': dt}
|
||||
# header_dict = Convert(data[next(_iter)][1:])
|
||||
header_dict.update(Convert(data[next(_iter)]))
|
||||
|
||||
# 读取波长等信息
|
||||
ls = []
|
||||
for i in range(int(header_dict['TotalSpectrometer'])):
|
||||
info = Convert(data[next(_iter)][1:]) # 略过第3行
|
||||
wvl = np.array(data[next(_iter)][1:])
|
||||
info['Wavelength'] = (wvl[wvl != '']).astype('float')
|
||||
info.update(header_dict)
|
||||
# info = my_dict
|
||||
ls.append(info)
|
||||
|
||||
# 读取光谱数据?????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
|
||||
data = data[next(_iter)+1:]
|
||||
data = pd.DataFrame(data)
|
||||
data[['spec', 'point']] = data[0].str.split('_', expand=True) #??
|
||||
data = data.set_index(['spec', 'point'])
|
||||
data = data.iloc[:, 3:]
|
||||
data.columns.name = 'Wavelength'
|
||||
data = data.astype('float')
|
||||
data = data.T.unstack().to_xarray()
|
||||
# _ ['spec'] = f
|
||||
data['Measures'] = dt
|
||||
data = data.expand_dims('Measures')
|
||||
return [data, ls] # data是光谱数据，ls是元信息
|
||||
|
||||
|
||||
def map_files(f):
|
||||
'''
|
||||
遍历文件
|
||||
'''
|
||||
dt = pd.to_datetime(f.split('\\')[-1][:-4], format='%Y_%m_%d_%H_%M_%S')
|
||||
data = read_files(f)
|
||||
data = data_parser(data, dt)
|
||||
# par_ls.append(par)
|
||||
return data # data是列表：第一个是光谱数据，第二个是元信息
|
||||
|
||||
# def map_files(i):
|
||||
# result = i * i
|
||||
# return result
|
||||
|
||||
|
||||
def processing(standard_sif, folder, out_file, pars, data, header, sky_p='P1', method='sfld'):
|
||||
'''
|
||||
提取算法调用
|
||||
'''
|
||||
# 参数解析为list
|
||||
pars = ast.literal_eval(pars)
|
||||
standard_sif = pd.read_csv(standard_sif, index_col=[0])
|
||||
point_ls = {}
|
||||
ls = []
|
||||
# 遍历光谱仪
|
||||
for spec in data['spec'].values:
|
||||
if "SIF" not in spec:
|
||||
continue
|
||||
else:
|
||||
#
|
||||
_ = data.sel(spec=spec) # ???????????????
|
||||
# 根据光谱仪提取对应的波长
|
||||
wvl = header['Wavelength'][spec]
|
||||
size = wvl.size
|
||||
# 使用波长的长度对数据进行截取（减少多光谱仪数据不一致导致的空值）
|
||||
_ = _.isel(Wavelength=xr.DataArray(
|
||||
np.arange(0, size), dims="Wavelength"))
|
||||
_['Wavelength'] = wvl
|
||||
# _ = _.where((_.Wavelength>731.3)&(_.Wavelength<782),drop=True)
|
||||
|
||||
sky = _.sel(point=sky_p, drop=True).rename('sky')
|
||||
for p in _.point:
|
||||
if p == sky_p:
|
||||
continue
|
||||
else:
|
||||
veg = _.sel(point=p, drop=True).rename('veg')
|
||||
input_each = xr.merge([sky, veg])
|
||||
_hf = intep(standard_sif, input_each.Wavelength.values) # 将标准sif插值匹配到数据的波长
|
||||
input_each['hf'] = (['Wavelength'], _hf)
|
||||
# 调用方法
|
||||
retr_method = sel_method(method)
|
||||
print('Running {} method on {} of spectrometer {}'.format(
|
||||
method, p.values, spec), flush=True)
|
||||
print('Processing ...', flush=True)
|
||||
sif = retr_method(input_each, *pars)[0]#---------------------------------------------------------------------------------------------
|
||||
point_ls.update({str(p.values): sif})
|
||||
_sif = pd.DataFrame(point_ls)
|
||||
point_ls = {}
|
||||
_sif.index = sky.Measures.values
|
||||
_sif.index.name = 'Measures'
|
||||
_sif = _sif.to_xarray()
|
||||
_sif = _sif.assign_coords(spec=spec)
|
||||
ls.append(_sif)
|
||||
sif = xr.merge(ls)
|
||||
sif.to_dataframe().to_csv(out_file)
|
||||
|
||||
|
||||
|
||||
|
||||
# # 需要输入的参数
|
||||
# # 标准SIF曲线
|
||||
# standard_sif = 'standard_sif.csv'
|
||||
# # 天空辐射对应的点
|
||||
# sky_p = 'P1'
|
||||
# # 存放文件夹
|
||||
# folder = '2021_12_21'
|
||||
# # 提取方法需要的输入参数，用括号括起来
|
||||
# pars = ([740,770],760)
|
||||
|
||||
# processing(standard_sif,folder,outfile,pars,sky_p='P1',method='svd')
|
||||
|
||||
def parse_args():
|
||||
parser = argparse.ArgumentParser(description="SIF retrieval process")
|
||||
parser.add_argument("standard_sif")
|
||||
parser.add_argument("folder")
|
||||
parser.add_argument("outfile")
|
||||
parser.add_argument("pars")
|
||||
parser.add_argument("sky_p", default='P1')
|
||||
parser.add_argument("method", default='sfld')
|
||||
args = parser.parse_args()
|
||||
return args
|
||||
|
||||
|
||||
def main():
|
||||
# 解析命令行参数
|
||||
inputs = parse_args()
|
||||
|
||||
# 读取数据
|
||||
cpus = os.cpu_count()
|
||||
pool = Pool(cpus)
|
||||
# print('tc-----------!', flush=True)
|
||||
files = glob(inputs.folder+'/*.csv')
|
||||
print('Total files: ', len(files), flush=True)
|
||||
|
||||
print('Reading files ...', flush=True)
|
||||
_ = pool.map(map_files, files)
|
||||
print('\tdone!', flush=True)
|
||||
|
||||
# 取出光谱数据
|
||||
# data = [d[0] for d in _]
|
||||
# for x in files:
|
||||
# map_files(files)
|
||||
|
||||
data = [d[0] for d in _]
|
||||
data = xr.concat(data, dim='Measures')
|
||||
# data.to_netcdf('C:\\fhz\\neibuwenjian\\tangchao\_Data\\temp.nc')
|
||||
|
||||
# 取出元数据
|
||||
header = _[0][1]
|
||||
_ = []
|
||||
header = pd.DataFrame(header).set_index('Model')
|
||||
# print(header, flush=True)
|
||||
|
||||
# sif算法处理
|
||||
processing(inputs.standard_sif, inputs.folder, inputs.outfile,
|
||||
inputs.pars, data, header, inputs.sky_p, inputs.method)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
1022
standard_sif.csv
Normal file
1022
standard_sif.csv
Normal file
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user