【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.mzph.cn/news/240641.shtml

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

相关文章

通过windows cng api 实现rsa非对称加密

参考&#xff1a; 1,使用 CNG 加密数据 - Win32 apps | Microsoft Learn 2,不记得了 &#xff08;下文通过cng api演示rsa加密&#xff0c;不做原理性介绍&#xff09; 相对于aes等对称加密算法&#xff0c;rsa加密算法不可逆性更强。非对称加密在通常情况下&#xff0c;使…

前端传输formDate格式的数据,后端不能用@RequestBody接收

写了个接口&#xff0c;跟前端对接&#xff0c;前端说怎么一直415的报错 我寻思不对啊&#xff0c;我swagger都请求成功了&#xff0c;后来发现前端一直是以formdata格式提交的数据&#xff0c;这样我其实是可以不加RequestBody的&#xff1b; 知识点&#xff1a; RequestBody…

类和对象

1 类定义&#xff1a; class ChecksumAccumulator {// class definition goes here } 你就能创建 ChecksumAccumulator 对象&#xff1a;new CheckSumAccumulator 注&#xff1a;1scala类中成员默认是public类型&#xff0c;若设为私有属性则必须加private关键字。在scala中是…

基于Springboot的留守儿童爱心网站(有报告)。Javaee项目,springboot项目。

演示视频&#xff1a; 基于Springboot的留守儿童爱心网站(有报告)。Javaee项目&#xff0c;springboot项目。 项目介绍&#xff1a; 采用M&#xff08;model&#xff09;V&#xff08;view&#xff09;C&#xff08;controller&#xff09;三层体系结构&#xff0c;通过Spring…

技术交底二维码的应用

二维码技术交底可以逐级落实、责任到人、有据可查、是目前最方便、实用的交底方式&#xff0c;下面我们讲解技术交底二维码的应用。 1、生成对应的技术交底二维码&#xff0c;将施工方案、技术资料、安全教育资料等内容上传到二维码里。打印出来现场粘贴&#xff0c;便于作业班…

java中线程相关的面试题

什么是线程安全&#xff0c;造成线程安全的本质是什么&#xff1f; 什么是线程安全呢&#xff1f; 咱们初步去理解话记住一句话就行&#xff1a;如果一个对象可以安全地被多个线程同时使用&#xff0c;那它就是线程安全的。 为什么并发编程会导致线程不安全&#xff1f; 可见…

用友U8CRM系统help2 任意文件读取漏洞复现

用友U8CRM系统的help2文件中接口存在任意文件读取漏洞&#xff0c;攻击者在未登录情况下即可进行漏洞利用。 1.1 漏洞级别 高危 1.2 快速检索 fofa语法&#xff1a; title"用友U8CRM"1.3 漏洞复现 该漏洞利用非常简单&#xff0c;只需构造get请求 访问该地址即可…

青少年CTF-qsnctf-A1-Misc-签到

题目环境&#xff1a; 题目难度&#xff1a;★题目描述&#xff1a;有没有可能&#xff0c;这个平台就是个题目&#xff1f; 一道杂项题 题目说的是这个平台就是题目 那么也就是说flag就在这个平台里面1.从高层次向低层次逐一排查 2.首先对平台首页进行排查进平台首页 第一种解…

HTML---浮动

提示&#xff1a;文章写完后&#xff0c;目录可以自动生成&#xff0c;如何生成可参考右边的帮助文档 文章目录 前言一、pandas是什么&#xff1f;二、使用步骤 1.引入库2.读入数据总结 一.常见的网页布局 二.标准文档流 标准文档流常见标签 标准文档流的组成 块级元素<div…

二叉搜索树 --- C++实现

目录 1.二叉搜索树的概念 2.二叉搜索树的操作 3. 二叉树的实现 4.二叉搜索树的应用 5. 二叉树的性能分析 6. 二叉树进阶练习题 1.二叉搜索树的概念 二叉搜索树又称二叉排序树&#xff0c;它或者是一棵空树&#xff0c;或者是具有以下性质的二叉树&#xff1a; 若它的左…

JDBC 知识点总结篇

JDBC 知识点总结篇 JDBC 接口 Java DataBase Connectivity Java数据库连接&#xff0c;由官方定义的一套操作所有关系型数据库的规则&#xff0c;即接口&#xff0c;各个数据库厂商实现该套接口 代码 // 本代码只提供一个样例&#xff0c;请根据自己实际情况修改代码 // 1.…

MyBatis笔记

Mybatis Mybatis介绍 什么是Mybatis? mybatis是支持普通SQL查询、存储过程和高级映射的优秀持久层框架。 Mybatis优点 几乎消除了JDBC代码和参数的手动设置消除结果集的检索使用XML或注解用于配置和原始映射&#xff0c;将接口和POJOs(实体类)映射成数据库中的记录。 My…

vue微乾坤子应用开发及ele组件开发时问题记录

一. 微乾坤 1. 新增page页面路由,pmi权限中心配置正常&#xff0c;跳转链接正确&#xff0c;但路由未找到403. 解决&#xff1a; 新增的配置是page类型&#xff0c;transformQianKunRoute方法转换微前端路由数据 时&#xff0c;过滤未兼容page型的路由&#xff0c; 解决 [menu,…

react中使用redux最简单最方便的方式,配合rematch简化操作,5分钟学会

react中使用状态管理的方式也很多&#xff0c;比如redux和mobx等&#xff0c;今天这一片就讲一下redux的入门到熟练使用&#xff0c;主要是要理解它redux的组成有哪些&#xff0c;到怎么创建&#xff0c;和组建中怎么使用三个问题。这里先放上官网文档&#xff0c;不理解的地方…

十一.约束(二)

约束 5.自增列:AUTO_INCREMENT5.1作用5.2关键字5.3特点和要求5.4如何指定自增约束5.5如何删除自增列5.6MySQL8.0新特性——自增变量的持久化 6.FOREIGN KEY 约束6.1作用6.2关键字6.3主表和从表/父表和子表6.4特点6.5添加外键约束6.6演示问题6.7约束等级6.8删除外键约束6.9开发场…

鸿蒙开发者工具安装及入门程序

下载工具DevEco Studio IDE 官网下载&#xff1a;HUAWEI DevEco Studio和SDK下载和升级 | HarmonyOS开发者 开发工具的安装 解压下载好的压缩包&#xff0c;一路无脑安装即可&#xff0c;安装完的使用方法类似于IDEA、WebStorm的使用&#xff0c;快捷键一致&#xff0c;默认黑…

【笔记】Spring的循环依赖

Spring的循环依赖 ObjectFactory:函数式接口&#xff0c;可以将lambda表达式作为参数放在方法的实参种&#xff0c;在方法执行的时候&#xff0c;并不会实际的调用当前lambda表达式&#xff0c;只有在调用getObject方法的时候才回去调用lambda表达式 为什么spring要用三级缓存…

常用的百兆网络变压器与RJ45网口的参考连接电路有哪些,主要注意事项在哪里呢?

Hqst华轩盛(石门盈盛)电子导读&#xff1a;一起来了解常用的百兆网络变压器与RJ45网口的参考连接电路有哪些&#xff0c;主要注意事项在哪里呢&#xff1f; 第一,常用的百兆网络变压器与RJ45网口的参考连接电路 常用百兆网络变压器与网口连接器分开为独立电子元件的分离式参考电…

TrustZone之与非安全虚拟化交互

到目前为止&#xff0c;我们在示例中忽略了非安全状态中可能存在的虚拟化程序。当存在虚拟化程序时&#xff0c;虚拟机与安全状态之间的许多通信将通过虚拟化程序进行。 例如&#xff0c;在虚拟化环境中&#xff0c;SMC用于访问固件功能和可信服务。固件功能包括诸如电源管理之…

Ubuntu 常用命令之 top 命令用法介绍

&#x1f4d1;Linux/Ubuntu 常用命令归类整理 top命令是Linux下常用的性能分析工具&#xff0c;可以实时动态地查看系统中各个进程的资源占用状况&#xff0c;类似于Windows的任务管理器。它可以显示系统总的和分区的CPU使用率、内存使用率、交换区使用率、系统负载、进程数、…