本次我们将分别使用两个流程完成对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
- Bit 31: Adjacency correction performed
| |
Bitmask for sur_refl_state_500m
|
- 0: Clear
- 1: Cloudy
- 2: Mixed
- 3: Not set, assumed clear
- 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
- Bit 11: Internal fire algorithm flag
- Bit 12: MOD35 snow/ice flag
- Bit 13: Pixel is adjacent to cloud
- Bit 14: BRDF correction performed data
- Bit 15: Internal snow mask
| |
数据属性:
这里再次提示大家,我们如果上传自己的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
对影像集合先去云后合成