【Earth Engine】协同Sentinel-1/2使用随机森林回归实现高分辨率相对财富（贫困）制图

目录

  • 1 简介与摘要
  • 2 思路
  • 3 效果预览
  • 4 代码思路
  • 5 完整代码
  • 6 后记

1 简介与摘要

最近在做一些课题，需要使用Sentinel-1/2进行机器学习制图。
然后想着总结一下相关数据和方法，就花半小时写了个代码。
然后再花半小时写下这篇博客记录一下。
因为基于多次拍脑袋而且花的时间很少，所以千万不要把这篇博客的实验流程当作完整的实验，
具体的科研实验还需要反复设计和多次试错！！！

这篇博客主要参考数据（相对贫困指数，RWI）来自这个GEE社区网站，有缺数据的可以直接在这上面找，要用的时候调用一下就行了。这是这个数据的一个交互地图预览。
工作完全是在GEE平台上写的，如上面所说，这个工作跟我的课题内容关系不大，纯粹是拍脑袋的需求然后拍脑袋写的代码。

2 思路

思路就是用Sentinel-1/2的一些波段作为自变量，对自变量在2400m进行采样（因为这个RWI数据的分辨率就是2400），
然后相对贫困指数RWI作为因变量，
最后在10m尺度使用随机森林回归进行相对贫困制图。

因为我没相关需求，求快，所以这里就简单随便弄弄。
因为是随便弄弄，我也不计算什么乱七八糟的指数了，就用VV和VH还有一些光学波段作为自变量。
因变量就直接用原数据的RWI。

研究区选的是坦桑尼亚的原首都达雷斯萨拉姆及其周边，选这个地方没啥特殊的原因，纯粹是点开那个交互地图映入眼帘的就是这个地方（好耳熟的名字，依稀记得当时所里好像有好多非洲老哥来自这个地方）。

所以总的来说，就是基于Sentinel-1/2做一个高分辨率的（10m）相对贫困地图。

这个博客主要还是展示如何用开源数据在GEE上快速实现机器学习制图，这个博客提供的思路太简单了，如果正经是搞科研论文或者做实验还涉及很多数据预处理步骤。
首先这个RWI数据就要咔咔一顿处理，在这个思路里直接在2400m采样肯定是行不通的。
所以这个博客主要还是提供一个参考，具体的一些流程还需要精心设计和反复试错。

3 效果预览

惯例，在代码前面先放效果预览。

我的兴趣区（roi）：
在这里插入图片描述

控制台：
在这里插入图片描述
上图的意思是区域内有2578个样本点，我拿80%（2078个）的去训练，20%（500个）去验证。
散点图是观测值和预测值的散点图。
RMSE就是均方根误差了。看着效果还可以哈。

绘制结果：
在这里插入图片描述

红色的是RWI高值（有钱），蓝色绿色的是RWI低值（贫困）。

绘制结果的细节：
在这里插入图片描述
在这里插入图片描述

为什么看起来这么稀碎呢，因为我拿World Settlement Footprint, WSF掩膜了一下，常规操作，具体可看代码。

4 代码思路

首先我把我要用的数据拉进来。其中RWI数据我按我的ROI筛选一下，不然太多了。代码如下：

// Importing Built-up map
var wsf2019 = ee.ImageCollection('projects/sat-io/open-datasets/WSF/WSF_2019');var wsf_01 = wsf2019.mosaic().eq(255);
var buildingup = wsf_01;// Importing RWI
var rwi = ee.FeatureCollection("projects/sat-io/open-datasets/facebook/relative_wealth_index").filterBounds(roi.geometry());print("sample size", rwi.size())

考虑到要调用Sentinel-1和2，所以也把这些数据写进来，然后也要写上一些预处理。代码如下：

// Importing S1 SAR images.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD').filterDate("2019-01-01", "2024-01-01").filterBounds(roi.geometry());
// print(sentinel1)
// Filter the S1 collection by metadata properties.
var vvVhIw = sentinel1// Filter to get images with VV and VH dual polarization..filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV')).filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))// Filter to get images collected in interferometric wide swath mode..filter(ee.Filter.eq('instrumentMode', 'IW'));
// Separate ascending and descending orbit images into distinct collections.
var vvVhIwAsc = vvVhIw.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));
var vvVhIwDesc = vvVhIw.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));// Importing S2 images.
// Cloud mask function.
function maskS2clouds(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));return image.updateMask(mask).divide(10000);
}var sentinel2 = ee.ImageCollection("COPERNICUS/S2_SR").filterDate("2019-01-01", "2024-01-01").filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)).filterBounds(roi.geometry()).map(maskS2clouds).mean();

再然后，把我要的波段挑出来，合并成一个image。简简单单搞一下子。代码如下：

// Indice.
var vv = vvVhIwAsc.merge(vvVhIwDesc).select('VV').mean();
var vh = vvVhIwAsc.merge(vvVhIwDesc).select('VH').mean();var blue = sentinel2.select('B2');
var green = sentinel2.select('B3');
var red = sentinel2.select('B4');
var red_edge1 = sentinel2.select('B5');
var red_edge2 = sentinel2.select('B6');
var red_edge3 = sentinel2.select('B7');
var nir = sentinel2.select('B8');
var red_edge4 = sentinel2.select('B8A');
var SWIR1 = sentinel2.select('B11');
var SWIR2 = sentinel2.select('B12');var image = ee.Image().addBands(vv).addBands(vh).addBands(blue).addBands(green).addBands(red).addBands(red_edge1).addBands(red_edge2).addBands(red_edge3).addBands(nir).addBands(red_edge4).addBands(SWIR1).addBands(SWIR2);var band_name = ee.List(['VV', 'VH', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12'])

然后，就是用RWI数据对我选的这些波段（自变量进行采样）。考虑到RWI是2400m分辨率，所以我就在2400米采样。采样完把我们样本规模打印到控制台看看。代码如下：

// Sampling.
var samples = rwi.randomColumn('random'); // set a random fieldvar training_samples = samples.filter(ee.Filter.lt('random', 0.8)); // 80% training data
var training_samples = image.reduceRegions({collection: training_samples, reducer: ee.Reducer.mean(), scale: 2400}); // samplingvar testing_samples = samples.filter(ee.Filter.gte('random', 0.8)); // 20% testing data
var testing_samples = image.reduceRegions({collection: testing_samples, reducer: ee.Reducer.mean(), scale: 2400}); // samplingprint("training sample size", training_samples.size())
print("testing sample size", testing_samples.size())

接下来就是要正式的训练模型绘图。用刚才采样的训练集训练一个回归模型，并且用这个模型进行绘图，得到一张结果的图像rwi_map，然后再用WSF掩膜一下住区。代码如下：

// Regressing.
var regressor = ee.Classifier.smileRandomForest(50) // number of estimators/trees.setOutputMode('REGRESSION') // regression algorithm.train(training_samples, // training samples"rwi", // regressing targetband_name); // regressor featuresvar rwi_map = image.unmask(0).clip(roi).classify(regressor).rename('rwi'); // regress
var rwi_map = rwi_map.updateMask(buildingup).clip(roi); // mask

建模画图完了，然后就是验证。用sampleRegions在画完的图上采样得到预测值，和一开始的观测值进行对比。验证包含散点图R2和RSME，都是必要常用的指标和表现形式。代码如下：

// Testing.
var testing_samples = rwi_map.rename('predicted').sampleRegions({collection: testing_samples,scale: 2400,geometries: true,projection: 'EPSG:4326'
});
var testing_samples = testing_samples.select(['rwi', 'predicted']);
// print(testing_samples)
// Scatter plot.
var chartTraining = ui.Chart.feature.byFeature({features: testing_samples, xProperty:'rwi', yProperties:['predicted']}).setChartType('ScatterChart').setOptions({title: 'Predicted vs Observed - Training data ',hAxis: {'title': 'observed'},vAxis: {'title': 'predicted'},pointSize: 3,trendlines: {0: {showR2: true,visibleInLegend: true},1: {showR2: true,visibleInLegend: true}}});
print(chartTraining);
// Calculating RMSE.
var observation = ee.Array(testing_samples.aggregate_array('rwi'));
var prediction = ee.Array(testing_samples.aggregate_array('predicted'));
var residuals = observation.subtract(prediction);
var rmse = residuals.pow(2).reduce('mean', [0]).sqrt();
print('RMSE', rmse);

最后是把刚才画的图，也就是绘图结果可视化。代码如下：

// Visualization.
Map.centerObject(roi, 8);
var palettes = ["#08525e", "#40d60e", "#b9e40e", "#f9c404", "e70000"]
Map.addLayer(rwi_map,{min: -0.5, max: 1.5, palette: palettes},'RWI');

5 完整代码

我的代码链接在这，可以直接使用。
完整代码如下（和链接中相同）：

// Importing Built-up map
var wsf2019 = ee.ImageCollection('projects/sat-io/open-datasets/WSF/WSF_2019');var wsf_01 = wsf2019.mosaic().eq(255);
var buildingup = wsf_01;// Importing RWI
var rwi = ee.FeatureCollection("projects/sat-io/open-datasets/facebook/relative_wealth_index").filterBounds(roi.geometry());print("sample size", rwi.size())// Importing S1 SAR images.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD').filterDate("2019-01-01", "2024-01-01").filterBounds(roi.geometry());
// print(sentinel1)
// Filter the S1 collection by metadata properties.
var vvVhIw = sentinel1// Filter to get images with VV and VH dual polarization..filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV')).filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))// Filter to get images collected in interferometric wide swath mode..filter(ee.Filter.eq('instrumentMode', 'IW'));
// Separate ascending and descending orbit images into distinct collections.
var vvVhIwAsc = vvVhIw.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));
var vvVhIwDesc = vvVhIw.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));// Importing S2 images.
// Cloud mask function.
function maskS2clouds(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));return image.updateMask(mask).divide(10000);
}var sentinel2 = ee.ImageCollection("COPERNICUS/S2_SR").filterDate("2019-01-01", "2024-01-01").filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)).filterBounds(roi.geometry()).map(maskS2clouds).mean();// Indice.
var vv = vvVhIwAsc.merge(vvVhIwDesc).select('VV').mean();
var vh = vvVhIwAsc.merge(vvVhIwDesc).select('VH').mean();var blue = sentinel2.select('B2');
var green = sentinel2.select('B3');
var red = sentinel2.select('B4');
var red_edge1 = sentinel2.select('B5');
var red_edge2 = sentinel2.select('B6');
var red_edge3 = sentinel2.select('B7');
var nir = sentinel2.select('B8');
var red_edge4 = sentinel2.select('B8A');
var SWIR1 = sentinel2.select('B11');
var SWIR2 = sentinel2.select('B12');var image = ee.Image().addBands(vv).addBands(vh).addBands(blue).addBands(green).addBands(red).addBands(red_edge1).addBands(red_edge2).addBands(red_edge3).addBands(nir).addBands(red_edge4).addBands(SWIR1).addBands(SWIR2);var band_name = ee.List(['VV', 'VH', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12'])// Sampling.
var samples = rwi.randomColumn('random'); // set a random fieldvar training_samples = samples.filter(ee.Filter.lt('random', 0.8)); // 80% training data
var training_samples = image.reduceRegions({collection: training_samples, reducer: ee.Reducer.mean(), scale: 2400}); // samplingvar testing_samples = samples.filter(ee.Filter.gte('random', 0.8)); // 20% testing data
var testing_samples = image.reduceRegions({collection: testing_samples, reducer: ee.Reducer.mean(), scale: 2400}); // samplingprint("training sample size", training_samples.size())
print("testing sample size", testing_samples.size())// Regressing.
var regressor = ee.Classifier.smileRandomForest(50) // number of estimators/trees.setOutputMode('REGRESSION') // regression algorithm.train(training_samples, // training samples"rwi", // regressing targetband_name); // regressor featuresvar rwi_map = image.unmask(0).clip(roi).classify(regressor).rename('rwi'); // regress
var rwi_map = rwi_map.updateMask(buildingup).clip(roi); // mask// Testing.
var testing_samples = rwi_map.rename('predicted').sampleRegions({collection: testing_samples,scale: 2400,geometries: true,projection: 'EPSG:4326'
});
var testing_samples = testing_samples.select(['rwi', 'predicted']);
// print(testing_samples)
// Scatter plot.
var chartTraining = ui.Chart.feature.byFeature({features: testing_samples, xProperty:'rwi', yProperties:['predicted']}).setChartType('ScatterChart').setOptions({title: 'Predicted vs Observed - Training data ',hAxis: {'title': 'observed'},vAxis: {'title': 'predicted'},pointSize: 3,trendlines: {0: {showR2: true,visibleInLegend: true},1: {showR2: true,visibleInLegend: true}}});
print(chartTraining);
// Calculating RMSE.
var observation = ee.Array(testing_samples.aggregate_array('rwi'));
var prediction = ee.Array(testing_samples.aggregate_array('predicted'));
var residuals = observation.subtract(prediction);
var rmse = residuals.pow(2).reduce('mean', [0]).sqrt();
print('RMSE', rmse);// Visualization.
Map.centerObject(roi, 8);
var palettes = ["#08525e", "#40d60e", "#b9e40e", "#f9c404", "e70000"]
Map.addLayer(rwi_map,{min: -0.5, max: 1.5, palette: palettes},'RWI');

6 后记

可能有些地方不太专业不太科学，希望诸位同行前辈不吝赐教或者有什么奇妙的想法可以和我共同探讨一下。可以在csdn私信我或者联系我邮箱（chinshuuichi@qq.com），不过还是希望大家邮箱联系我，csdn私信很糟糕而且我上csdn也很随缘。
如果对你有帮助，还望支持一下~点击此处施舍或扫下图的码。
-----------------------分割线（以下是乞讨内容）-----------------------
在这里插入图片描述

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.hqwc.cn/news/294461.html

如若内容造成侵权/违法违规/事实不符，请联系编程知识网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

浅谈Guava Cache的参数使用

浅谈Guava Cache的参数使用

CacheLoader 用于数据加载方式比较固定且统一的场景&#xff0c;在缓存容器创建的时候就需要指定此具体的加载逻辑。通常开发中使用时我们需要继承CacheLoader类或写一个匿名实现类实现其load方法和reload方法 load方法 当执行get操作没有命中缓存或者判断缓存已经超出expir…
阅读更多...
NiNNet

NiNNet

目录 一、网络介绍 1、全连接层存在的问题 2、NiN的解决方案(NiN块) 3、NiN架构 4、总结 二、代码实现 1、定义NiN卷积块 2、NiN模型 3、训练模型 一、网络介绍 NiN&#xff08;Network in Network&#xff09;是一种用于图像识别任务的卷积神经网络模型。它由谷歌研究…
阅读更多...
Linux学习小结

Linux学习小结

目录结构 tree -L 1 / # /root #root用户的家目录 /home #存储普通用户家目录 lostfound #这个目录平时是空的&#xff0c;存储系统非正常关机而留下“无家可归”的文件 /usr #系统文件&#xff0c;相当于C:\Windows /usr/local #软件安装的目录&#xff0c;相当于C:\Progra…
阅读更多...
跨境电商独立站深度分析演示网站

跨境电商独立站深度分析演示网站

对于跨境电商卖家来说&#xff0c;多平台、多站点的布局是非常重要的战略。这样做可以规避”鸡蛋放在同一个篮子里”的风险也能够追求更高的销售额和利润。同时&#xff0c;市场的变化也带来了新的发展机会&#xff0c;因此很多出海企业都希望抓住独立站的新机遇&#xff0c;抢…
阅读更多...
【华为数据之道学习笔记】6-4 打造数据供应的“三个1”

【华为数据之道学习笔记】6-4 打造数据供应的“三个1”

数据服务改变了传统的数据集成方式&#xff0c;所有数据都通过服务对外提供&#xff0c;用户不再直接集成数据&#xff0c;而是通过服务获取。因此&#xff0c;数据服务应该拉动数据供应链条的各个节点&#xff0c;以方便用户能准确地获取数据为重要目标。 数据供应到消费的完整…
阅读更多...
【Linux笔记】文件和目录操作

【Linux笔记】文件和目录操作

&#x1f34e;个人博客&#xff1a;个人主页 &#x1f3c6;个人专栏&#xff1a;Linux学习 ⛳️ 功不唐捐&#xff0c;玉汝于成 目录 前言 命令 ls (List): pwd (Print Working Directory): cp (Copy): mv (Move): rm (Remove): 结语 我的其他博客 前言 学习Linux命令…
阅读更多...
开发知识点-HTML/JavaScript

开发知识点-HTML/JavaScript

HTML/JavaScript xlinksvgviewBoxuse基础预热与语法基础知识js 如何运行页面适用js 及输出 面向对象抽奖功能 json 支持 字符串转数组数组转字符串数组元素删除长度0位添加一个元素// 表示在下标为1处添加一项tttarray.splice(1,0,ttt)//[123,ttt,456]// 数组是否包含某个元素a…
阅读更多...
PIC单片机项目（8）——基于PIC16F877A的温度光照检测装置的protues仿真

PIC单片机项目（8）——基于PIC16F877A的温度光照检测装置的protues仿真

1.功能设计 使用PIC16F877A单片机&#xff0c;进行温度检测、光照检测。温度使用的是DS18B20&#xff0c;光照检测直接利用的AD转换。 光照太暗就开灯&#xff0c;温度太高就开风扇。温度阈值和光照阈值都实时显示在LCD1602屏幕上面。 完成了protues仿真。文件里面包含代码和仿…
阅读更多...
blender径向渐变材质-着色编辑器

blender径向渐变材质-着色编辑器

要点&#xff1a; 1、用纹理坐标中的物体输出连接映射中的矢量输入 2、物体选择一个空坐标&#xff0c;将空坐标延z轴上移一段距离 3、空坐标的大小要缩放到和要添加材质的物体大小保持一致
阅读更多...
微前端样式隔离、sessionStorage、localStorage隔离

微前端样式隔离、sessionStorage、localStorage隔离

1、样式隔离 前端样式不隔离&#xff0c;会产生样式冲突的问题&#xff0c;这个点在qiankun也存在 子应用1修改一个样式 button {background: red&#xff01;important&#xff1b; }其它应用也会受到影响 qiankun的css隔离方案&#xff08;shadow dom&#xff09; shadow …
阅读更多...
华为 1+X 网络系统运维与建设中级实操模拟题

华为 1+X 网络系统运维与建设中级实操模拟题

目 实验拓扑 配置中的注意事项&#xff1a;&#xff08;针对新手&#xff09; 实验目的 实验要求 实验步骤 一、搭建实验拓扑 二、配置主机名称 三、配置链路聚合 四、VLAN 配置 五、配置 RSTP 协议 六、配置 IP 地址 七、配置 VRRP 协议。 八、配置 OSPF 协议 九…
阅读更多...
C语言中常用的sscanf函数

C语言中常用的sscanf函数

文章目录 1. 接受全部参数&#xff1a;2、分辨数字和字符3. 数字和字符一起会默认是字符4. 同时接收多个变量5. 指定长度的集合操作6. 排除部分字符 sscanf()定义于头文件stdio.h。sscanf()会将参数str的字符串根据参数format字符串来转换并格式化数据。格式转换形式请参考scan…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023