PIE-engine 教程 ——MODIS影像去云教程(山西省为例)

本次我们将分别使用两个流程完成对MODIS影像去除云，第一个就是先去云然后再合成，第二个方式是先合成后去云，我们通常情况下一般都是先去云再合成。

本文使用的数据：

MOD09A1产品来自于MODIS Terra星，数据为经过了大气校正（如气体，气溶胶和瑞利散射）的地表光谱反射率的估计值。数据包含7个反射率波段，及2个质量波段和4个观测波段，分辨率为500米。
依据云层和太阳天顶情况，从8天中选择一个最优值作为像素值。当存在多个最优值时，选取蓝波段（band 3）最小的像素值作为8天合成的值。

名称	单位	类型	分辨率(m)	波长(nm)	比例因子	值域范围	无效值	描述信息
sur_refl_b01	--	int16	500	620~670	0.0001	-100~16000	-28672	波段1的表面反射率
sur_refl_b02	--	int16	500	841~876	0.0001	-100~16000	-28672	波段2的表面反射率
sur_refl_b03	--	int16	500	459~479	0.0001	-100~16000	-28672	波段3的表面反射率
sur_refl_b04	--	int16	500	545~565	0.0001	-100~16000	-28672	波段4的表面反射率
sur_refl_b05	--	int16	500	1230~1250	0.0001	-100~16000	-28672	波段5的表面反射率
sur_refl_b06	--	int16	500	1628~1652	0.0001	-100~16000	-28672	波段6的表面反射率
sur_refl_b07	--	int16	500	2105~2155	0.0001	-100~16000	-28672	波段7的表面反射率
sur_refl_qc_500m	--	int16	--	--	--	--		表面反射率500m波段质量控制标志
sur_refl_szen	Degrees	int16	--	--	0.01	0~18000	0	太阳天顶角
sur_refl_vzen	Degrees	int16	--	--	0.01	0~18000	0	视角天顶角
sur_refl_raz	Degrees	int16	--	--	0.01	-18000~18000	0	相对方位角
sur_refl_state_500m	--	uint16	--	--	--	--	65535	表面反射率500m状态标志
sur_refl_day_of_year	--	uint16	--	--	--	1~366	65535	儒略日(在儒略周期内以连续的日数计算时间的计时法)

Bitmask for sur_refl_qc_500m

Bits 0-1: MODLAND QA bits 0: Corrected product produced at ideal quality - all bands 1: Corrected product produced at less than ideal quality - some or all bands 2: Corrected product not produced due to cloud effects - all bands 3: Corrected product not produced for other reasons - some or all bands, may be fill value (11) [Note that a value of (11) overrides a value of (01)] Bit 2-5: Band 1 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bits 6-9: Band 2 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bits 10-13: Band 3 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bits 14-17: Band 4 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bits 18-21: Band 5 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bits 22-25: Band 6 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bits 26-29: Band 7 data quality, four bit range 0: Highest quality 7: Noisy detector 8: Dead detector, data interpolated in L1B 9: Solar zenith ≥ 86 degrees 10: Solar zenith ≥ 85 and < 86 degrees 11: Missing input 12: Internal constant used in place of climatological data for at least one atmospheric constant 13: Correction out of bounds, pixel constrained to extreme allowable value 14: L1B data faulty 15: Not processed due to deep ocean or clouds Bit 30: Atmospheric correction performed 0: No 1: Yes Bit 31: Adjacency correction performed 0: No 1: Yes

Bitmask for sur_refl_state_500m

Bits 0-1: Cloud state 0: Clear 1: Cloudy 2: Mixed 3: Not set, assumed clear Bit 2: Cloud shadow 0: No 1: Yes Bits 3-5: Land/water flag 0: Shallow ocean 1: Land 2: Ocean coastlines and lake shorelines 3: Shallow inland water 4: Ephemeral water 5: Deep inland water 6: Continental/moderate ocean 7: Deep ocean Bits 6-7: Aerosol quantity 0: Climatology 1: Low 2: Average 3: High Bits 8-9: Cirrus detected 0: None 1: Small 2: Average 3: High Bit 10: Internal cloud algorithm flag 0: No cloud 1: Cloud Bit 11: Internal fire algorithm flag 0: No fire 1: Fire Bit 12: MOD35 snow/ice flag 0: No 1: Yes Bit 13: Pixel is adjacent to cloud 0: No 1: Yes Bit 14: BRDF correction performed data 0: No 1: Yes Bit 15: Internal snow mask 0: No snow 1: Snow

数据属性：

date

string

影像日期

这里再次提示大家，我们如果上传自己的featureCollection想要进行筛选和裁剪影像的时候我们需要完成对矢量集合的处理，.first().geometry()来完成准换。另外，这里我们镶嵌使用的函数：

mosaic()

将影像集合融合成为一张影像Image，融合规则保留是这个影像集合中最新的有效像素值。

方法参数：

- imageCollection(ImageCollection)

ImageCollection实例。

返回值：Image

代码：

var featureCollection0 = pie.FeatureCollection('user/xg346049806/shanxibianjie').first().geometry();

//这个是和GEE不同的一点，就是需要先获取一个提取bit的函数
function bitwiseExtract(value, fromBit, toBit) {
  if (toBit === undefined) toBit = fromBit;
  var maskSize = 1 + toBit - fromBit; //位数
  var mask = pie.Number(1 << maskSize).subtract(1); //3=11，即挑选位
  return value.rightShift(fromBit).bitwiseAnd(mask);
}
//如果在GEE中就可以直接使用这个函数就可以了
function cloudfree_mod09a1(image) {
  //提取QA波段
  var qa = image.select("sur_refl_state_500m");
  //提取云比特位
  var cloudState = bitwiseExtract(qa, 0, 1);
  //提取云阴影比特位
  var cloudShadowState = bitwiseExtract(qa, 2, 2);
  //提取卷云比特位
  var cirrusState = bitwiseExtract(qa, 8, 9);
  //叠加得到无云、无云阴影且无卷云的区域
  var mask = cloudState
    .eq(0) // Clear
    .and(cloudShadowState.eq(0)) // No cloud shadow
    .and(cirrusState.eq(0)); // No cirrus
  return image.updateMask(mask);
}

//筛选影像
var imgcol_m = pie
  .ImageCollection("USGS/MOD09A1/006")
  .filterDate("2019-01-01", "2019-02-01")
  .filterBounds(featureCollection0)
  .select([
    "sur_refl_b01",
    "sur_refl_b04",
    "sur_refl_b03",
    "sur_refl_state_500m",
  ]);
print("imgcol_m", imgcol_m);

//（1）对影像集合先合成后去云
var raw_img_m = imgcol_m.mosaic().clip(featureCollection0);
print("raw_img_m", raw_img_m);
Map.addLayer(
  raw_img_m,
  {
    min: 0,
    max: 3500,
    bands: ["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"],
  },
  "raw_img_m"
);
Map.centerObject(geometry0, 6);
var result = cloudfree_mod09a1(raw_img_m);
Map.addLayer(
  result,
  {
    min: 0,
    max: 3500,
    bands: ["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"],
  },
  "result"
);
//（2）对影像集合先去云后合成
var cf_img_m = imgcol_m.map(cloudfree_mod09a1).mosaic().clip(featureCollection0);
Map.addLayer(
  cf_img_m,
  {
    min: 0,
    max: 3500,
    bands: ["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"],
  },
  "cf_img_m"
);

对影像集合先合成后去云 

 result

 对影像集合先去云后合成