当前位置：首页 > news >正文

基于 GEE 利用 RF 回归模型实现空间降尺度

news 来源：原创 2025/5/12 20:45:43

1 前言

2 完整代码

3 运行结果

1 前言

本篇讲述在GEE上基于回归模型降尺度，也就是要复现以下论文，该论文发表在J-Star期刊上。

“Ebrahimy H, Aghighi H, Azadbakht M, et al. Downscaling MODIS land surface temperature product using an adaptive random forest regression method and Google Earth Engine for a 19-years spatiotemporal trend analysis over Iran[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 2103-2112.”

这篇文章的整体流程图如下：

其核心部分就是蓝色框内的内容，整体上和PCA降尺度大同小异，主要涉及到的操作就是准备数据、数据重采样、函数映射规则构建和降尺度结果。其中，文章说的climate zone可以忽略不计，因为这是一个外部文件。也就是说，类似于气候分区，这一步是预先就有了的。文章只是把这一点放大了，说是adaptive和climate-aware的，其实整个流程并不涉及，所以整体上并没有多大效果。在复现的时候不考虑这个。

那么，基于随机森林回归（RF regression）模型的降尺度和基于PCA的降尺度的区别就在于模型的构建，其实和分类也差不多，只不过分类最后对应的是离散量，而回归对应的是连续变量。此外，和PCA降尺度对比可以发现，PCA最后是线性模型，也就是构建了LST与PCA分量之间的关系，而求解的过程中是逐元素求解的，所以之前很多人不理解。而基于RF归回的降尺度，从操作的角度来说更简单，因为它只需要抽样就行了，然后建立回归，最后将模型映射到不同尺度数据上就可以了。

选择一个研究区进行测试，最后的结果如下。可以发现，经过校正后的LST，分辨率更高了，同时也保持了原有的空间分布特征。

2 完整代码

var roi = table;
Map.centerObject(roi);
Map.addLayer(roi, {'color':'grey'}, 'studyarea');var day_start = '2018-05-01';
var day_end = '2018-06-01';// MODIS daily LST product (MOD11A1)
var MOD_LST = ee.ImageCollection('MODIS/006/MOD11A1').filter(ee.Filter.date(day_start, day_end)).select(['LST_Day_1km'],['LST']).map(function(img) {return img.multiply(0.02).clip(roi)});
var MOD_LST = MOD_LST.select('LST').median();
// print("MOD_LST",MOD_LST);var LSTVisPara = {min: 290.0,max: 310.0,palette: ['040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6','0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef','3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f','fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d','ff0000', 'de0101', 'c21301', 'a71001', '911003'],
};
Map.addLayer(MOD_LST, LSTVisPara,'MOD_LST',false);// MODIS daily surface reflectance product (MOD09GQ)
var MOD_SR = ee.ImageCollection('MODIS/006/MOD09GQ').filter(ee.Filter.date(day_start, day_end)).select(['sur_refl_b01','sur_refl_b02'],['red','NIR']).map(function(img){var ndvi = img.normalizedDifference(['NIR','red']).rename("NDVI")return img.addBands(ndvi).clip(roi);}).median();var falseColorVis = {min: -100.0,max: 8000.0,bands: ['NIR', 'NIR', 'red'],
};
Map.addLayer(MOD_SR, falseColorVis, 'MOD_SR',false);// MODIS yearly land cover product (MCD12Q1) 
var MOD_LULC = ee.ImageCollection('MODIS/006/MCD12Q1').select(['LC_Type1'],['LULC']).filter(ee.Filter.date('2018-01-01','2019-01-01')).map(function(img) {return img.clip(roi)}).mean();var igbpLandCoverVis = {min: 0.0,max: 17.0,palette: ['05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159','dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c','69fff8', 'f9ffa4', '1c0dff'],
};
Map.addLayer(MOD_LULC, igbpLandCoverVis, 'MOD_LULC',false);// the Shuttle Radar Topography Mission (SRTM)
var SRTM = ee.Image('CGIAR/SRTM90_V4').select('elevation').clip(roi);
// Map.addLayer(SRTM, {min: 0, max: 300}, 'SRTM',false);var img_Original = MOD_LST.addBands(MOD_SR).addBands(MOD_LULC).addBands(SRTM);
print("img_Original",img_Original);var img_240 = reSampleFun(img_Original,240);
var img_960 = reSampleFun(img_Original,960);
Map.addLayer(img_240.select("LST"), LSTVisPara,'img_240',false);
Map.addLayer(img_960.select("LST"), LSTVisPara,'img_960',false);var randomPoint = ee.FeatureCollection.randomPoints({region:roi,points:5000
});
Map.addLayer(randomPoint, {'color':'black'}, 'randomPoint');var SamplePointCol =  img_960.sampleRegions({collection: randomPoint,   // properties: ['agbValue','植被类'],scale: 100,tileScale: 2,geometries:true
});
print("SamplePointCol",SamplePointCol.limit(100));var classifierOfChoice = ee.Classifier.smileRandomForest(500).setOutputMode('REGRESSION');
var propertyToPredictAsString = 'LST';
var bandFeature = ee.List(img_960.bandNames()).remove('LST');
var trainedModel = classifierOfChoice.train({features: SamplePointCol,inputProperties: bandFeature,classProperty: propertyToPredictAsString
});    
var LST960 = img_960.select(bandFeature)                   .classify(trainedModel)                   .rename('LST960').clip(roi);
Map.addLayer(LST960.select("LST960"), LSTVisPara,'LST960',false); var LSTError = img_960.select("LST").subtract(LST960.select("LST960"));
var LSTError_240 = reSampleFun(LSTError,240);
var LST240 = img_240.select(bandFeature)                   .classify(trainedModel)                   .rename('LST240').clip(roi);
Map.addLayer(LST240.select("LST240"), LSTVisPara,'LST240',false); var LST240_corrected  = LST240.select("LST240").add(LSTError_240).rename("LST240_corrected");
Map.addLayer(LST240_corrected.select("LST240_corrected"), LSTVisPara,'LST240_corrected',false); // resample image
function reSampleFun(image,scale){var dst_crs = image.select('blue').projection().crs();var reSampleImg = image // .resample('bilinear') // bilinear bicubic.reproject({crs:'EPSG:4326', scale:scale});return reSampleImg;
}

3 运行结果

最后，进行一些简单的说明：

这里只是实现了最基本的操作，很多细节部分并没有涉及到，其中最主要的就是如何抽选出最具有代表性的样本数据来构建模型，这一点很关键，所以大家如果想使用的话，还是需要根据自己的研究进行适当调整。
理论上来说，这里已经实现了降尺度，但是这种方式还需要进一步完善。

1 前言

2 完整代码

3 运行结果

相关文章：