线性插值
- 1. 写在前面
- 2. MODIS数据插值
- 3. Landsat数据插值
- 3.1 参数修改以适应其他类型的遥感数据
- 3.2 Landsat数据汇总
- 3.3 Sentinel卫星介绍
1. 写在前面
之前我写了使用年内均值或者中值来填补数据控制的方法，这种方法较为简单，不够精确。因此，我在这里结合所学知识提供一种线性插值的方法。线性插值假设在两个已知数据点之间的数据变化是均匀的，通过直线连接这两个点，并据此估算中间点的值。这种方法在数据缺失较少且缺失点两侧都有有效数据时效果较好。
本研究以成都市成华区作为感兴趣区，选择MODIS和Landsat8数据进行线性插值。
研究区：
2. MODIS数据插值
MOD09Q1产品只有红光（sur_refl_b01）和近红外（sur_refl_b02）两个波段的地表反射率数据。但是空间分辨率为250m，相对较高。
这种方法是从相关公众号学到的，代码如下：
// GS（2024）0650号========================================================================================
var china_provinces = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_provinces");
var china_city = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_city");
var china_county = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_county");
var sichuan = china_provinces.filter(ee.Filter.eq('name','四川省'))
var chengdu = china_city.filter(ee.Filter.eq('name','成都市'))
var chenhuaqu = china_county.filter(ee.Filter.eq('name','成华区'))
Map.addLayer(sichuan.style({fillColor:'00000000',color:'red'}),{},"四川省", 0) //ff0000
Map.addLayer(chengdu.style({fillColor:'00000000',color:'blue'}),{},"成都市", 1) //ffff00
Map.addLayer(chenhuaqu.style({color:"black"}),{},"成华区")
Map.centerObject(chenhuaqu, 11.7);
var colorizedVis = {
min: -1.0,
max: 1.0,
palette: [
'0000ff','0000ee','0000dd','0000cc','0000bb','0000aa','000099',
'000088','000077','000066','000055','000044','ce7e45', 'df923d',
'f1b555', 'fcd163', '99b718','74a901', '66a000', '207401', '012e01',
'262600', '0e0e00'
]
};
//参数设置;==============================================================================================
var mod_sr = ee.ImageCollection("MODIS/061/MOD09Q1"); //2000-02-18T00:00:00Z–2024-09-05T00:00:00Z
var startDate = ee.Date("2023-06-01");
var endDate = ee.Date("2023-09-01");
var roi = chenhuaqu
function replace_mask(img, newimg, nodata) {
nodata = nodata || 0;
var mask = img.mask();
img = img.unmask(nodata);
img = img.where(mask.not(), newimg);
img = img.updateMask(img.neq(nodata));
return img;
}
/** Interpolation not considering weights */
function addTimeBand(img) {
/** make sure mask is consistent */
var mask = img.mask();
var time = img.metadata('system:time_start').rename("time").mask(mask);
return img.addBands(time);
}
var getQABits = function(image, start, end) {
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += 1 << i;
}
// Return a single band image of the extracted QA bits
return image.bitwiseAnd(pattern).rightShift(start);
};
function linearInterp(imgcol, frame, nodata){
frame = frame || 32;
nodata = nodata || 0;
var time = 'system:time_start';
imgcol = imgcol.map(addTimeBand);
var maxDiff = ee.Filter.maxDifference(frame * (1000*60*60*24), time, null, time);
var cond = {leftField:time, rightField:time};
// Images after, sorted in descending order (so closest is last).
var f1 = ee.Filter.and(maxDiff, ee.Filter.lessThanOrEquals(cond));
var c1 = ee.Join.saveAll({matchesKey:'after', ordering:time, ascending:false})
.apply(imgcol, imgcol, f1);
// Images before, sorted in ascending order (so closest is last).
var f2 = ee.Filter.and(maxDiff, ee.Filter.greaterThanOrEquals(cond));
var c2 = ee.Join.saveAll({matchesKey:'before', ordering:time, ascending:true})
.apply(c1, imgcol, f2);
var interpolated = ee.ImageCollection(c2.map(function(img) {
img = ee.Image(img);
var before = ee.ImageCollection.fromImages(ee.List(img.get('before'))).mosaic();
var after = ee.ImageCollection.fromImages(ee.List(img.get('after'))).mosaic();
img = img.set('before', null).set('after', null);
// constrain after or before no NA values, confirm linear Interp having result
before = replace_mask(before, after, nodata);
after = replace_mask(after , before, nodata);
// Compute the ratio between the image times.
var x1 = before.select('time').double();
var x2 = after.select('time').double();
var now = ee.Image.constant(img.date().millis()).double();
var ratio = now.subtract(x1).divide(x2.subtract(x1)); // this is zero anywhere x1 = x2
// Compute the interpolated image.
before = before.select(0); //remove time band now;
after = after.select(0);
img = img.select(0);
var interp = after.subtract(before).multiply(ratio).add(before).toFloat();
var qc = img.mask().not().rename('qc');
interp = replace_mask(img, interp, nodata).rename('MOD_NDVI_INTER');
return interp.addBands(qc).copyProperties(img, img.propertyNames());
}));
return interpolated;
}
function main(){
var mod_ndvi = mod_sr.filterDate(startDate,endDate).filterBounds(roi)
.map(function(image){return image.clip(roi)})
.map(function(image) {
var ndvi = image.normalizedDifference(['sur_refl_b02', 'sur_refl_b01']); //calculate Landsat NDVI
var Q = image.select('State');
var masknum = getQABits(Q, 0, 2);
var mask = masknum.eq(ee.Image(0)).or(masknum.eq(ee.Image(3)));
var ndvi_masked = ndvi.updateMask(mask); //generate MODIS cloud mask from state band (0:clear;1: cloud/shadow/mixed)
return ndvi_masked.rename(['MOD_NDVI']).copyProperties(image, image.propertyNames())});
var nodata = -9999;
var frame = 4*4;
var mod_ndvi_interp = linearInterp(mod_ndvi,frame,nodata);
mod_ndvi = mod_ndvi.toList(100)
mod_ndvi_interp = mod_ndvi_interp.select('MOD_NDVI_INTER').toList(100);
print('mod_ndvi_interp',mod_ndvi_interp)
Map.addLayer(ee.Image(mod_ndvi.get(5)),colorizedVis,'mod_ndvi')
Map.addLayer(ee.Image(mod_ndvi_interp.get(5)),colorizedVis,'mod_ndvi_interp')
}
main();
结果对比：
3. Landsat数据插值
在这里我在MODIS数据插值的基础上进行了修改和调整，并且添加了注释以便读者阅读。我这里使用了Landsat8数据为例，并且我会提供如何修改这份代码的方法，使得读者能够实现对其他类型的遥感数据进行插值，如sentinel卫星数据。
代码如下：
//*********************************************************************************************************
//*********************************************************************************************************
// GS（2024）0650号========================================================================================
var china_provinces = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_provinces");
var china_city = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_city");
var china_county = ee.FeatureCollection("projects/ee-tilmacatanla/assets/boundry/china_county");
var sichuan = china_provinces.filter(ee.Filter.eq('name','四川省'))
var chengdu = china_city.filter(ee.Filter.eq('name','成都市'))
var chenhuaqu = china_county.filter(ee.Filter.eq('name','成华区'))
Map.addLayer(sichuan.style({fillColor:'00000000',color:'red'}),{},"四川省", 0) //ff0000
Map.addLayer(chengdu.style({fillColor:'00000000',color:'blue'}),{},"成都市", 1) //ffff00
Map.addLayer(chenhuaqu.style({color:"black"}),{},"成华区")
Map.centerObject(chenhuaqu, 11.7);
var colorizedVis = {
min: -1.0,
max: 1.0,
palette: [
'0000ff','0000ee','0000dd','0000cc','0000bb','0000aa','000099',
'000088','000077','000066','000055','000044','ce7e45', 'df923d',
'f1b555', 'fcd163', '99b718','74a901', '66a000', '207401', '012e01',
'262600', '0e0e00'
]
};
var colorizedVis1 = {
min: -1.0,
max: 1.0,
palette: ['#d7191c','#fdae61','#ffffbf','#a6d96a','#1a9641']
};
// 参数设置================================================================================================
var startDate = ee.Date("2023-03-01");
var endDate = ee.Date("2023-08-01");
//函数定义=================================================================================
// Applies scaling factors.
function applyScaleFactors(image) {
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
return image.addBands(opticalBands, null, true)
.addBands(thermalBands, null, true);
}
function rmCloudNew(image) {
var cloudShadowBitMask = (1 << 4);
var cloudsBitMask = (1 << 3);
var qa = image.select('QA_PIXEL');
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return image.updateMask(mask)
.copyProperties(image)
.copyProperties(image, ["system:time_start"]);
}
// Landsat8===================================================================================
var mod_sr = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterBounds(chenhuaqu)
.filter(ee.Filter.calendarRange(2014,2023,'year'))
.filter(ee.Filter.calendarRange(1,12,'month'))
.map(applyScaleFactors)
.map(rmCloudNew);
// 将一幅图像的掩膜区域用另一幅图像的对应区域替换，同时确保所有 nodata 值的像素都被掩膜====================
function replace_mask(img, newimg, nodata) {
nodata = nodata || 0;
var mask = img.mask();
img = img.unmask(nodata);
img = img.where(mask.not(), newimg);
img = img.updateMask(img.neq(nodata));
return img;
}
/** Interpolation not considering weights */
function addTimeBand(img) {
/** make sure mask is consistent */
var mask = img.mask();
var time = img.metadata('system:time_start').rename("time").mask(mask);
return img.addBands(time);
}
var getQABits = function(image, start, end) {
var pattern = 0;
for (var i = start; i <= end; i++) {
pattern += 1 << i;
}
// Return a single band image of the extracted QA bits
return image.bitwiseAnd(pattern).rightShift(start);
};
function linearInterp(imgcol, Interp_range, nodata){
Interp_range = Interp_range || 30;
nodata = nodata || 0;
var time = 'system:time_start';
imgcol = imgcol.map(addTimeBand);
print("imgcol:", imgcol)
//选择图像集合中时间差异不超过 Interp_range 天（转换为毫秒数）的图像
var maxDiff = ee.Filter.maxDifference(Interp_range * (1000*60*60*24), time, null, time);//Interp_range 通常以天为单位, 这里将天数转换为毫秒数
var cond = {leftField:time, rightField:time};
// Images after, sorted in descending order (so closest is last).
var f1 = ee.Filter.and(maxDiff, ee.Filter.lessThanOrEquals(cond)); // 比较图像自身的时间字段，用于确保图像是按照时间顺序排列
var c1 = ee.Join.saveAll({
matchesKey:'after', //寻找时间上在当前图像之后的图像
ordering:time, //定义匹配图像时的排序字段，这里是time，即按照图像获取的时间进行排序
ascending:false //定义了排序的方向，false表示按降序排序，意味着最接近当前图像时间的图像将被排在最后
}).apply(imgcol, imgcol, f1);
print("c1:", c1)
// Images before, sorted in ascending order (so closest is last).
var f2 = ee.Filter.and(maxDiff, ee.Filter.greaterThanOrEquals(cond));
var c2 = ee.Join.saveAll({
matchesKey:'before', //寻找时间上在当前图像之前的图像
ordering:time, //定义匹配图像时的排序字段，这里是time，即按照图像获取的时间进行排序
ascending:true //定义了排序的方向，true表示按升序排序
}).apply(c1, imgcol, f2);
print("c2:", c2)
var interpolated = ee.ImageCollection(
c2.map(function(img) {
img = ee.Image(img);
// 获取roi图像缺失数据点前后最接近的图像，以便根据这些图像来估计缺失数据点的值
var before = ee.ImageCollection.fromImages(ee.List(img.get('before'))).mosaic();// mosaic获取的是最后一个有值的像素，然后生成一张影像；
var after = ee.ImageCollection.fromImages(ee.List(img.get('after'))).mosaic();
// 属性初始化
img = img.set('before', null).set('after', null);
// constrain after or before no NA values, confirm linear Interp having result
// 更新 before 图像，使其掩膜区域被 after 图像的对应区域替换，同时确保所有 nodata 值的像素都被掩膜
before = replace_mask(before, after, nodata);
after = replace_mask(after , before, nodata);
// Compute the ratio between the image times.
var x1 = before.select('time').double();
var x2 = after.select('time').double();
var now = ee.Image.constant(img.date().millis()).double();
var ratio = now.subtract(x1).divide(x2.subtract(x1)); // this is zero anywhere x1 = x2
// Compute the interpolated image.
before = before.select(0); //remove time band now; 选择 before 图像的第一个波段
after = after.select(0);
img = img.select(0);
var interp = after.subtract(before).multiply(ratio).add(before).toFloat();
var qc = img.mask().not().rename('qc');
interp = replace_mask(img, interp, nodata).rename('Landsat_NDVI_INTER');
return interp.addBands(qc).copyProperties(img, img.propertyNames());
}));
return interpolated;
}
function main(){
var Landsat_NDVI = mod_sr.filterDate(startDate,endDate).filterBounds(chenhuaqu)
.map(function(image){return image.clip(chenhuaqu)})
.map(function(image) {
var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']); //calculate Landsat NDVI
var Q = image.select('QA_PIXEL'); // Pixel quality
// 从原始图像中筛选出无云或云状态未设置的像素
// 提取第8-11位(不可行)===============================================================
// Bit 0: Fill; Bit 1: Dilated Cloud; Bit 2: Cirrus (high confidence);
// Bit 3: Cloud; Bit 4: Cloud Shadow; Bit 5: Snow; Bit 6: Clear
// Bits 8-9:Cloud Confidence：云层置信度，两位用来表示云层存在的可能性，从无到高
// Bits 10-11: Cloud Shadow Confidence:云阴影置信度，两位用来表示云阴影存在的可能性，从无到高
var masknum = getQABits(Q, 0, 4);
var mask = masknum.eq(ee.Image(0)).or(masknum.eq(ee.Image(64))); // 0001000000 = 64~[2*(7-1)=64]
var ndvi_masked = ndvi.updateMask(mask); //generate Landsat cloud mask from QA_PIXEL band (0:clear;1: cloud/shadow/mixed)
return ndvi_masked.rename(['Landsat_NDVI']).copyProperties(image, image.propertyNames())
});
//=====================================================================================
var imgTiles = ee.FeatureCollection(Landsat_NDVI.map(function(img){
var tmpFootprint = ee.Image(img).geometry();
var WRS_PATH = ee.Number(ee.Image(img).get("WRS_PATH"));
var WRS_ROW = ee.Number(ee.Image(img).get("WRS_ROW"));
return ee.Feature(tmpFootprint,null).set("WRS_PATH",WRS_PATH).set("WRS_ROW",WRS_ROW);
})).distinct(["WRS_PATH",'WRS_ROW']);
var styling = {color:'red',fillColor:'00000000'};
print("imgTiles size",imgTiles);
Map.addLayer(imgTiles.style(styling), {}, 'imgTiles');
// 计算每个图像的研究区像元数并添加为属性=============================================
var withPixelCount = Landsat_NDVI.map(function(image) {
// 假设我们的NDVI波段名为'NDVI'
var pixelCount = image.reduceRegion({
reducer: ee.Reducer.count(),
geometry: chenhuaqu,
scale: 900, // 根据实际需要调整分辨率
maxPixels: 1e9
});
// 使用波段名称作为键来获取计数
return image.set('pixelCount', pixelCount.get('Landsat_NDVI'));
});
// 根据像元数对图像集合进行排序
var sortedCollection = withPixelCount.sort('pixelCount', true);
// 展示排序后的图像集合。
var count = sortedCollection.size().getInfo();
for (var i = 0; i < count; i++) {
var image = ee.Image(sortedCollection.toList(count).get(i));
var dateString = image.get('system:index').getInfo();
var pixelCount = image.get('pixelCount').getInfo();
Map.addLayer(image.clip(chenhuaqu), colorizedVis1, 'Image: ' + dateString + ', Pixels: ' + pixelCount, 0);
}
var nodata = -9999;
var Interp_range = 60; //*******************************此处可以修改为你自己的天数***************************
var Landsat_NDVI_interp = linearInterp(Landsat_NDVI,Interp_range,nodata);
Landsat_NDVI = Landsat_NDVI.toList(100)
var numberOfImagesInList = Landsat_NDVI.size();
print('Number of images in list:', numberOfImagesInList);
// 计算每个图像的研究区像元数并添加为属性=============================================
var withPixelCount_Inter = Landsat_NDVI_interp.map(function(image) {
// 假设我们的NDVI波段名为'NDVI'
var pixelCount = image.reduceRegion({
reducer: ee.Reducer.count(),
geometry: chenhuaqu,
scale: 900, // 根据实际需要调整分辨率
maxPixels: 1e9
});
// 使用波段名称作为键来获取计数
return image.set('pixelCount', pixelCount.get('Landsat_NDVI_INTER'));
});
// 根据像元数对图像集合进行排序
var sortedCollection_Inter = withPixelCount_Inter.sort('pixelCount', true);
// 展示排序后的图像集合。
var count_Inter = sortedCollection_Inter.size().getInfo();
for (var j = 0; j < count_Inter; j++) {
var image_Inter = ee.Image(sortedCollection_Inter.toList(count_Inter).get(j));
var dateString_Inter = image_Inter.get('system:index').getInfo();
var pixelCount_Inter = image_Inter.get('pixelCount').getInfo();
Map.addLayer(image_Inter.select('Landsat_NDVI_INTER').clip(chenhuaqu), colorizedVis, 'Inter_Image: ' + dateString_Inter + ', Pixels: ' + pixelCount_Inter, 0);
Export.image.toDrive({
image: image_Inter,
description: 'Landsat_NDVI_Interp_' + dateString_Inter,
folder: 'NDVI_Export',
scale: 30,
region: chenhuaqu,
maxPixels: 1e9,
});
}
Landsat_NDVI_interp = Landsat_NDVI_interp.select('Landsat_NDVI_INTER').toList(100);
print('Landsat_NDVI_interp',Landsat_NDVI_interp)
}
main();
结果对比：
由此可见插值前后发生了巨大变化。
c1结果：
c2结果：
在原有的基础上，我增加了显示遥感数据条带号的函数，这里可以结合我之前写个一篇博客进行学习：GEE26：遥感数据可用数据源计算及条带号制作
结果如下：
3.1 参数修改以适应其他类型的遥感数据
这里我将介绍如何修改以上代码参数来增加代码的可用性。首先，我们需要对数据源进行修改：
var mod_sr = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterBounds(chenhuaqu)
.filter(ee.Filter.calendarRange(2014,2023,'year'))
.filter(ee.Filter.calendarRange(1,12,'month'))
.map(applyScaleFactors)
.map(rmCloudNew);
此外，你需要修改function main()函数。首先，这里需要调整用于计算NDVI的波段以及质量评价波段，Landsat使用的是SR_B5和SR_B4，并且使用QA_PIXEL波段评价数据质量：
function main(){
var Landsat_NDVI = mod_sr.filterDate(startDate,endDate).filterBounds(chenhuaqu)
.map(function(image){return image.clip(chenhuaqu)})
.map(function(image) {
var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']); //calculate Landsat NDVI
var Q = image.select('QA_PIXEL'); // Pixel quality
// 从原始图像中筛选出无云或云状态未设置的像素
// 提取第8-11位(不可行)===============================================================
// Bit 0: Fill; Bit 1: Dilated Cloud; Bit 2: Cirrus (high confidence);
// Bit 3: Cloud; Bit 4: Cloud Shadow; Bit 5: Snow; Bit 6: Clear
// Bits 8-9:Cloud Confidence：云层置信度，两位用来表示云层存在的可能性，从无到高
// Bits 10-11: Cloud Shadow Confidence:云阴影置信度，两位用来表示云阴影存在的可能性，从无到高
var masknum = getQABits(Q, 0, 4);
var mask = masknum.eq(ee.Image(0)).or(masknum.eq(ee.Image(64))); // 0001000000 = 64~[2*(7-1)=64]
var ndvi_masked = ndvi.updateMask(mask); //generate Landsat cloud mask from QA_PIXEL band (0:clear;1: cloud/shadow/mixed)
return ndvi_masked.rename(['Landsat_NDVI']).copyProperties(image, image.propertyNames())
});
QA_PIXEL介绍：
这里我使用了第7位，即Bit6来构建无云天气的掩膜数据，即0001000000 = 64。
再来看上面的MOD09Q1数据，使用的是：
var masknum = getQABits(Q, 0, 2);
var mask = masknum.eq(ee.Image(0)).or(masknum.eq(ee.Image(3)));
即00000011=3，而其对应的质量评估字段名称为state：
3.2 Landsat数据汇总
最后，我汇中了所有的Landsat数据，读者可以选择满足自己条件的数据进行以上操作：
3.3 Sentinel卫星介绍
- 哨兵－1卫星是全天时、全天候雷达成像任务，用于陆地和海洋观测，首颗哨兵－1A卫星已于2014年4月3日发射。
- 哨兵－2卫星是多光谱高分辨率成像任务，用于陆地监测，可提供植被、土壤和水覆盖、内陆水路及海岸区域等图像，还可用于紧急救援服务。
- 哨兵－3卫星携带多种有效载荷，用于高精度测量海面地形、海面和地表温度、海洋水色和土壤特性，还将支持海洋预报系统及环境与气候监测。
- 哨兵－4载荷专用于大气化学成分监测，将搭载在第三代气象卫星－S（MTG－S）上。
- 哨兵－5载荷用于监测大气环境，将搭载在欧洲第二代“气象业务”（MetOp）卫星上。
- 哨兵－5P卫星用于减小欧洲“环境卫星”（Envisat）和哨兵－5载荷之间的数据缺口。
- 哨兵－6卫星是贾森－3（Jason－3）海洋卫星的后续任务，将携带雷达高度计，用于测量全球海面高度，主要用于海洋科学和气候研究。
数据对比：
