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

news2024/11/24 15:39:03

目录

  • 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）：
在这里插入图片描述

控制台：
在这里插入图片描述

上图的意思是区域内有1943个样本点，我拿90%（1751个）的去训练，10%（192个）去验证。

绘制结果：

在这里插入图片描述

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

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

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

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'));
// Calculate temporal means.
// var vhIwAscMean = vvVhIwAsc.select('VH').mean();
// var vhIwDescMean = vvVhIwDesc.select('VH').mean();
// var vv = vvVhIwAsc.merge(vvVhIwDesc).select('VV').mean();
// var vh = vvVhIwAsc.merge(vvVhIwDesc).select('VH').mean();


// 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 field

var training_samples = samples.filter(ee.Filter.lt('random', 0.9)); // 90% training data
var training_samples = image.reduceRegions({collection: training_samples, reducer: ee.Reducer.mean(), scale: 2400}); // sampling


var testing_samples = samples.filter(ee.Filter.gte('random', 0.9)); // 10% testing data
var testing_samples = image.reduceRegions({collection: testing_samples, reducer: ee.Reducer.mean(), scale: 2400}); // sampling


print("training sample size", training_samples.size())
print("testing sample size", testing_samples.size())

接下来就是要正式的训练模型绘图。用刚才采样的训练集训练一个回归模型，并且用这个模型进行绘图，而且还可视化出来。代码如下：

// Regressing.
var regressor = ee.Classifier
            .smileRandomForest(20) // number of estimators/trees
            .setOutputMode('REGRESSION') // regression algorithm
            .train(training_samples, // training samples
                    "rwi", // regressing target
                    band_name); // regressor features


var rwi_map = image.unmask(0).clip(roi).classify(regressor).rename('rwi'); // regress
var rwi_map = rwi_map.updateMask(buildingup).clip(roi); // mask


// Visualization.
var palettes = ["#08525e", "#40d60e", "#b9e40e", "#f9c404", "e70000"]
Map.addLayer(rwi_map,
              {min: -0.2, max: 1.2, 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'));
// Calculate temporal means.
// var vhIwAscMean = vvVhIwAsc.select('VH').mean();
// var vhIwDescMean = vvVhIwDesc.select('VH').mean();
// var vv = vvVhIwAsc.merge(vvVhIwDesc).select('VV').mean();
// var vh = vvVhIwAsc.merge(vvVhIwDesc).select('VH').mean();


// 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 field

var training_samples = samples.filter(ee.Filter.lt('random', 0.9)); // 90% training data
var training_samples = image.reduceRegions({collection: training_samples, reducer: ee.Reducer.mean(), scale: 2400}); // sampling


var testing_samples = samples.filter(ee.Filter.gte('random', 0.9)); // 10% testing data
var testing_samples = image.reduceRegions({collection: testing_samples, reducer: ee.Reducer.mean(), scale: 2400}); // sampling


print("training sample size", training_samples.size())
print("testing sample size", testing_samples.size())


// Regressing.
var regressor = ee.Classifier
            .smileRandomForest(20) // number of estimators/trees
            .setOutputMode('REGRESSION') // regression algorithm
            .train(training_samples, // training samples
                    "rwi", // regressing target
                    band_name); // regressor features


var rwi_map = image.unmask(0).clip(roi).classify(regressor).rename('rwi'); // regress
var rwi_map = rwi_map.updateMask(buildingup).clip(roi); // mask


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

6 后记

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

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

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

二叉搜索树（AVL树，红黑树）+封装

二叉搜索树（AVL树，红黑树）+封装

就像学习其他的东西一样&#xff0c;首先我们要知道二叉搜索树的作用和定义是什么&#xff01; 首先顾名思义&#xff0c;二叉搜索树肯定是被用来为搜索服务的数据结构。 并且它的搜索效率可以达到logN,也就是一百万的数据也只用查找几十次&#xff08;AVL树可以控制在20次左…
阅读更多...
日常工作中常用的抓包工具都有哪些呢？

日常工作中常用的抓包工具都有哪些呢？

&#x1f4e2;专注于分享软件测试干货内容&#xff0c;欢迎点赞 &#x1f44d; 收藏 ⭐留言 &#x1f4dd; 如有错误敬请指正&#xff01;&#x1f4e2;交流讨论&#xff1a;欢迎加入我们一起学习&#xff01;&#x1f4e2;资源分享&#xff1a;耗时200小时精选的「软件测试」资…
阅读更多...
电子合同的分类有哪些？

电子合同的分类有哪些？

1、从电子合同订立的具体方式的角度&#xff0c;可分为利用电子数据交换订立的合同和利用电子邮件订立的合同&#xff1b; 2、从电子合同标的物的属性的角度&#xff0c;可分为网络服务合同、软件授权合同、需要物流配送的合同等&#xff1b; 3、从电子合同当事人的性质的角度…
阅读更多...
使用Gitee中的CI/CD来完成代码的自动部署与发布（使用内网穿透把本地电脑当作服务器使用）

使用Gitee中的CI/CD来完成代码的自动部署与发布（使用内网穿透把本地电脑当作服务器使用）

&#x1f4da;目录 &#x1f4da;简介:⚙️ 所需工具&#xff1a;&#x1f4a8;内网穿透配置&#x1f4ad;工具介绍✨命令安装&#x1f38a;配置Cpolar&#x1f573;️关闭防火墙&#x1f95b;防火墙端口放行规则&#xff08;关闭防火墙可以忽略&#xff09;&#x1f36c;小章总…
阅读更多...
【python】作用域与闭包 || global与nonlocal

【python】作用域与闭包 || global与nonlocal

python作用域 其他语言的作用域&#xff1a;块级、函数、类、模块、包等由小到大的级别但是python没有块级&#xff08;if语句块、for语句块&#xff09;&#xff0c;所以if中定义的变量&#xff0c;相当于普通语句 >>> if True: # if语句块没有作用域x …
阅读更多...
华为云之ECS云产品快速入门

华为云之ECS云产品快速入门

华为云之ECS云产品快速入门 一、ECS云服务器介绍二、本次实践目标三、创建虚拟私有云VPC1.虚拟私有云VPC介绍2.进入虚拟私有云VPC管理页面3.创建虚拟私有云4.查看创建的VPC 四、创建弹性云服务器ECS——Linux1.进入ECS购买界面2.创建弹性云服务器(Linux)——基础配置步骤3.创建…
阅读更多...
如何使用 templ 在 Go 中编写 HTML 用户界面？

如何使用 templ 在 Go 中编写 HTML 用户界面？

关注公众号【爱发白日梦的后端】分享技术干货、读书笔记、开源项目、实战经验、高效开发工具等&#xff0c;您的关注将是我的更新动力&#xff01; 简介 templ 是一个在 Go 中编写 HTML 用户界面的语言。使用 templ&#xff0c;我们可以创建可呈现 HTML 片段的组件&#xff0c…
阅读更多...
基于改进YOLOv7的绝缘子缺陷检测算法

基于改进YOLOv7的绝缘子缺陷检测算法

摘要 现有的检测方法面临着巨大的挑战&#xff0c;在识别绝缘子的微小缺陷时&#xff0c;针对输电线路图像与复杂的背景。为保证输电线路的安全运行&#xff0c;提出一种改进的YOLOv 7模型&#xff0c;以提高检测结果。 首先&#xff0c;基于K-means对绝缘子数据集的目标盒进…
阅读更多...
Redis“垃圾”过期死键管理与优化

Redis“垃圾”过期死键管理与优化

【作者】付磊 Redis死键的定义不尽相同&#xff0c;通常有两种&#xff1a; 写到Redis里后&#xff0c;由于过期时间过长或者压根没有过期时间&#xff0c;加之长期不访问&#xff0c;这类key可以被称为死键。 明明已经过了过期时间&#xff0c;但还占用Redis内存&#xff08…
阅读更多...
利用tf-idf对特征进行提取

利用tf-idf对特征进行提取

TF-IDF是一种文本特征提取的方法&#xff0c;用于评估一个词在一组文档中的重要性。 一、代码 from sklearn.feature_extraction.text import TfidfVectorizer import numpy as npdef print_tfidf_words(documents):"""打印TF-IDF矩阵中每个文档中非零值对应…
阅读更多...
Nacos-服务发现与配置管理v1.0

Nacos-服务发现与配置管理v1.0

Nacos - 服务发现和配置管理 教学目标 1&#xff09;能够理解微服务架构的特点 2&#xff09;能够理解服务发现的流程 3&#xff09;能够说出Nacos的功能 4&#xff09;掌握Nacos的安装方法 5&#xff09;掌握RESTful服务发现开发方法 6&#xff09;掌握Dubbo服务发现开…
阅读更多...
AIGC绘画关键词 - 神兽类（一）

AIGC绘画关键词 - 神兽类（一）

Unity3D特效百例案例项目实战源码Android-Unity实战问题汇总游戏脚本-辅助自动化Android控件全解手册再战Android系列Scratch编程案例软考全系列Unity3D学习专栏蓝桥系列ChatGPT和AIGC &#x1f449;关于作者 专注于Android/Unity和各种游戏开发技巧&#xff0c;以及各种资源分…
阅读更多...
JNI学习(二)

JNI学习(二)

静态注册 接着上篇博客学习 JNI函数 JNIEXPORT void JNICALL Java_com_example_jnidemo_TextDemo_setText(JNIEnv *env, jobject this, jstring string){ __android_log_print(ANDROID_LOG_ERROR, "test", "invoke set from C\n");char* str (char*)(*e…
阅读更多...
每日一题——链表的回文结构

每日一题——链表的回文结构

链表的回文结构 1. 题目描述 对于一个链表&#xff0c;请设计一个时间复杂度为O(n),额外空间复杂度为O(1)的算法&#xff0c;判断其是否为回文结构。 给定一个链表的头指针A&#xff0c;请返回一个bool值&#xff0c;代表其是否为回文结构。保证链表长度小于等于900。 测试…
阅读更多...
SpringBoot集成swagger-ui

SpringBoot集成swagger-ui

1.引入依赖&#xff1a; <!--swagger--><dependency><groupId>io.springfox</groupId><artifactId>springfox-swagger2</artifactId><version>2.7.0</version></dependency><dependency><groupId>io.sprin…
阅读更多...
回归预测 | MATLAB实现GWO-DHKELM基于灰狼算法优化深度混合核极限学习机的数据回归预测 （多指标，多图）

回归预测 | MATLAB实现GWO-DHKELM基于灰狼算法优化深度混合核极限学习机的数据回归预测 （多指标，多图）

回归预测 | MATLAB实现GWO-DHKELM基于灰狼算法优化深度混合核极限学习机的数据回归预测 &#xff08;多指标&#xff0c;多图&#xff09; 目录 回归预测 | MATLAB实现GWO-DHKELM基于灰狼算法优化深度混合核极限学习机的数据回归预测 &#xff08;多指标&#xff0c;多图&#…
阅读更多...
查找算法——二分查找

查找算法——二分查找

笔记&#xff1a;二分查找算法 | 数据结构与算法 系列教程&#xff08;笔记&#xff09; 题目描述 请对一个 有序数组 进行二分查找 {1,8, 10, 89, 1000, 1234}&#xff0c;输入一个数看看该数组是否存在此数&#xff0c;并且求出下 标&#xff0c;如果没有就提示「没有这个数…
阅读更多...
前端ICON库

前端ICON库

前端ICON库 1.mingcute mingcute 2.lordicon lordicon 3.字节iconpark&#xff08;推荐&#xff09; 字节iconpark 4.iconbuddy iconbuddy.app/ 5.商标寻找youicons 免费下载数百万个徽标以获得设计灵感 | YouIcons.com 还有一堆工具
阅读更多...
一键转换，将HTML智能转换为PDF，轻松解决文档转换需求

一键转换，将HTML智能转换为PDF，轻松解决文档转换需求

在数字时代&#xff0c;HTML网页是我们获取信息的主要来源之一。然而&#xff0c;有时候我们可能需要将网页内容以PDF格式保存&#xff0c;以便于离线阅读、打印或分享。这时&#xff0c;将HTML转换为PDF就变得尤为重要。 首先&#xff0c;我们要进入首助编辑高手主页面&#x…
阅读更多...
华为数通方向HCIP-DataCom H12-831题库(多选题：181-200）

华为数通方向HCIP-DataCom H12-831题库(多选题：181-200）

第181题 如图所示，R1、R2、R3、R4都部署为SPF区域0，链路的cost值如图中标识。R1、R2R3、R4的Loopback0通告入OSPF。R1、R2、R3与R4使用Loopback0作为连接接口,建立BGP对等体关系，其中R4为RR设备，R1、R2、R3是R4的客户端。当R4的直连地址172.20,1,4/32通告入BGP后，以下关R…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023