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 bands1: Corrected product produced at less than ideal quality - some or all bands2: Corrected product not produced due to cloud effects - all bands3: 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 quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 6-9: Band 2 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 10-13: Band 3 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 14-17: Band 4 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 18-21: Band 5 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 22-25: Band 6 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 26-29: Band 7 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBit 30: Atmospheric correction performed 0: No1: YesBit 31: Adjacency correction performed 0: No1: Yes

Bitmask for sur_refl_qc_500m

Bits 0-1: MODLAND QA bits 0: Corrected product produced at ideal quality - all bands1: Corrected product produced at less than ideal quality - some or all bands2: Corrected product not produced due to cloud effects - all bands3: 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 quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 6-9: Band 2 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 10-13: Band 3 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 14-17: Band 4 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 18-21: Band 5 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 22-25: Band 6 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBits 26-29: Band 7 data quality, four bit range 0: Highest quality7: Noisy detector8: Dead detector, data interpolated in L1B9: Solar zenith ≥ 86 degrees10: Solar zenith ≥ 85 and < 86 degrees11: Missing input12: Internal constant used in place of climatological data for at least one atmospheric constant13: Correction out of bounds, pixel constrained to extreme allowable value14: L1B data faulty15: Not processed due to deep ocean or cloudsBit 30: Atmospheric correction performed 0: No1: YesBit 31: Adjacency correction performed 0: No1: Yes

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: Clear1: Cloudy2: Mixed3: Not set, assumed clearBit 2: Cloud shadow 0: No1: YesBits 3-5: Land/water flag 0: Shallow ocean1: Land2: Ocean coastlines and lake shorelines3: Shallow inland water4: Ephemeral water5: Deep inland water6: Continental/moderate ocean7: Deep oceanBits 6-7: Aerosol quantity 0: Climatology1: Low2: Average3: HighBits 8-9: Cirrus detected 0: None1: Small2: Average3: HighBit 10: Internal cloud algorithm flag 0: No cloud1: CloudBit 11: Internal fire algorithm flag 0: No fire1: FireBit 12: MOD35 snow/ice flag 0: No1: YesBit 13: Pixel is adjacent to cloud 0: No1: YesBit 14: BRDF correction performed data 0: No1: YesBit 15: Internal snow mask 0: No snow1: Snow

Bitmask for sur_refl_state_500m

Bits 0-1: Cloud state 0: Clear1: Cloudy2: Mixed3: Not set, assumed clearBit 2: Cloud shadow 0: No1: YesBits 3-5: Land/water flag 0: Shallow ocean1: Land2: Ocean coastlines and lake shorelines3: Shallow inland water4: Ephemeral water5: Deep inland water6: Continental/moderate ocean7: Deep oceanBits 6-7: Aerosol quantity 0: Climatology1: Low2: Average3: HighBits 8-9: Cirrus detected 0: None1: Small2: Average3: HighBit 10: Internal cloud algorithm flag 0: No cloud1: CloudBit 11: Internal fire algorithm flag 0: No fire1: FireBit 12: MOD35 snow/ice flag 0: No1: YesBit 13: Pixel is adjacent to cloud 0: No1: YesBit 14: BRDF correction performed data 0: No1: YesBit 15: Internal snow mask 0: No snow1: Snow

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

代码：

代码语言：javascript

复制

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

对影像集合先去云后合成