1、质量控制：小于0的辐亮度值设置为很小的正值；

2、3fld和sfld：限制找OA位置的窗口到759-761nm；
3、3fld和sfld：限制找OA位置的数据的时间 → 早上10点到下午3点；

.gitignore

@@ -1,6 +1,8 @@
 # 唐超添加
 /.idea
-
+/problem
+/sample_data
+/paper
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

sif_methods.py

@@ -6,11 +6,16 @@ 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]
+    sky_spec = data.sky.to_pandas().between_time('10:00:00', '15:00:00')
+    veg_spec = data.veg.to_pandas().between_time('10:00:00', '15:00:00')
+
+    sky_idxmin = sky_spec.idxmin(axis=1)
+    veg_idxmin = veg_spec.idxmin(axis=1)
+
+    wvl_inside_band_l = veg_idxmin.mean()
+    wvl_inside_band_e = sky_idxmin.mean()
+
+    return[wvl_inside_band_l, wvl_inside_band_e]
 
 def cal_outside_values_mean(data,outer):
     '''
@@ -38,7 +43,9 @@ def sfld(data,wl_range,outer):
     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)
+
+    data2 = data.where((data.Wavelength > outer[1]) & (data.Wavelength < outer[1] + 3), drop=True)
+    [wvl_inside_band_l,wvl_inside_band_e]=cal_inside_bands_ave(data2)
     for i in range(0,nmeas_):
         _data = data.isel(Measures=i)
         veg_out,sky_out,_ = cal_outside_values_mean(_data,outer)
@@ -67,7 +74,9 @@ def fld3(data,wl_range,outer_left,outer_right):
     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)
+
+    data2 = data.where((data.Wavelength > outer_left[1]) & (data.Wavelength < outer_right[0]), drop=True)
+    [wvl_inside_band_l,wvl_inside_band_e]=cal_inside_bands_ave(data2)  # 应该局限到759-761
     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)
@@ -237,6 +246,10 @@ def doas(data,wl_range,band=760):
         print(data.Measures[i].values)
         sky_spec = data.sky[i].values
         veg_spec = data.veg[i].values
+
+        veg_spec[np.where(veg_spec < 0)] = 0
+        sky_spec[np.where(sky_spec < 0)] = 0
+
         _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_retrieval.py

@@ -65,7 +65,7 @@ def read_files(f):
     '''
     with open(f) as csvfile:
         data = list(csv.reader(csvfile))  # 读csv并转化为list
-#         data = np.array(data)
+        data = [[element for element in row if element != ''] for row in data]  # 去掉每行的空字符
         data = np.array(data, dtype=object)  # 转化为np数组
         data = [np.array(d) for d in data]
         csvfile.close()
@@ -123,11 +123,19 @@ 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是列表：第一个是光谱数据，第二个是元信息
+    try:
+        dt = pd.to_datetime(f.split('\\')[-1][:-4], format='%Y_%m_%d_%H_%M_%S')
+        data = read_files(f)  # 每个list都是一个np数组
+        data = data_parser(data, dt)
+    except Exception as e:
+        print("Reading error: %s" % f)
+    else:  # 如果打开没有出错
+        # par_ls.append(par)
+        return data    # data是列表：第一个是光谱数据，第二个是元信息
+    finally:
+        pass
+
+
 
 # def map_files(i):
 #     result = i * i
@@ -160,11 +168,21 @@ def processing(standard_sif, folder, out_file, pars, data, header, sky_p='P1', m
         #     _ = _.where((_.Wavelength>731.3)&(_.Wavelength<782),drop=True)
 
             sky = _.sel(point=sky_p, drop=True).rename('sky')
+            # 质量控制
+            tmp = sky.values
+            tmp[tmp < 0] = 0.0000000000000001
+            sky.values = tmp
+
             for p in _.point:
                 if p == sky_p:
                     continue
                 else:
                     veg = _.sel(point=p, drop=True).rename('veg')
+                    # 质量控制
+                    tmp = veg.values
+                    tmp[tmp < 0] = 0.0000000000000001
+                    veg.values = tmp
+
                     input_each = xr.merge([sky, veg])
                     _hf = intep(standard_sif, input_each.Wavelength.values)  # 将标准sif插值匹配到数据的波长
                     input_each['hf'] = (['Wavelength'], _hf)
@@ -281,4 +299,5 @@ def main():
 
 
 if __name__ == '__main__':
+    print('Version: 2.3.6')
     main()