Google Earth Engine（GEE）——Sentinel-1 和 2 数据的融合，水稻范围识别和水稻种植季节区分地图绘制—马来西亚为例
最近发现了一个巨牛的人工智能学习网站，点击跳转到网站：前言 – 床长人工智能教程
这次给大家推荐一篇文章，关于水稻识别：使用 GEE 平台中的无监督分类整合 Sentinel-1 和 2 时间序列数据来实现的。这项研究的结果将为绘制稻田及其生长阶段的地图提供一个有效的框架，这对于解决粮食安全问题和减少水稻种植产生的甲烷气体排放具有重要价值。
摘要：
水稻是世界一半以上人口的主要作物，但缺乏概述水稻产区及其生长阶段的高分辨率地图。大多数遥感研究绘制了水稻的范围：然而，在热带地区，水稻全年种植，种植日期和种植频率不同。因此，绘制水稻生长阶段比仅绘制范围更，有用。本研究通过开发一种基于物候学的方法解决了这一挑战。假设是Sentinel-1和2时间序列数据的无监督分类（k均值聚类）可以识别稻田和生长阶段，因为（1）Sentinel-1 VH反向散射可以识别移栽过程中是否存在洪水；（2）稻田生，长阶段（营养、生殖、和成熟阶段）直到收获点可以通过归一化差异植被指数（NDVI）时间序列来识别。使用所提出的方法，本研究绘制了马来西亚半鸟（131，598公里）的稻田范围和种植日历2）在Google地球引攀（GEE）平台上。使用k-means聚类对2019年1月至2020年12月的Sentinel-1和2月度时间序列数据进行分类，以识别具有相似，物候模式的区域。这种方法产生了10米分辨率的稻田范围、强度和种植日历图。使用来自谷歌地球的高分辨率街景图像进行验证表明，预测地图的总体准确率为95.95%，kappa系数为0.92，此外，预测的农作物日历与当地政府的粮仓数据吻合良好。结果表明，所提出的基于物候学的方法具有成本效益，并且可以准确地绘制大面积稻田和生长阶段的地图。
关键词：稻田：种植模式：物候学：哨兵1：哨兵2：谷歌地球引擎：机器学习：深度学习
为了提取由 K-means 聚类识别的具有代表性的 VH 极化和 NDVI 聚类单元，每个聚类随机生成总共 1000 个样本。随后，每个集群的每月时间序列的空间中值被用作代表性的 VH 极化和 NDVI 时间序列剖面。
Sentinel-1数据（VH极化值）对检测移栽期更敏感；另一方面，Sentinel-2 数据（NDVI 值）更详细地显示了生长和成熟阶段。局部最小值和峰值之间的差异约为四个月，代表一个生长季节。休耕期 VH 极化和 NDVI 值下降。
还有一个APP：
malaysiarice
文中的技术流程：
本文需要用到的函数：
toBands()
将一个集合转换为一个单一的多波段图像，包含集合中每个图像的所有波段。输出的带子是通过在现有的带子名称前加上它来自的图像ID来命名的（例如：'image1_band1'）。
参数。
this:集合（ImageCollection）。
输入的集合。
返回。图像
sample(region, scale, projection, factor, numPixels, seed, dropNulls, tileScale, geometries)
对图像的像素进行采样，并将其作为一个特征集合返回。每个特征在输入图像中的每一带都有一个属性。请注意，默认行为是放弃与屏蔽像素相交的特征，这将导致空值属性（见dropNulls参数）。
参数。
this:image（图像）。
要采样的图像。
region (Geometry, default: null):
要取样的区域。如果没有指定，则使用图像的整个面积。
scale（Float，默认为空）。
要取样的投影的名义比例，以米为单位。
projection（投影，默认为空）。
要采样的投影。如果没有指定，则使用图像的第一个波段的投影。如果与scale同时指定，则按指定的比例重新调整。
factor (Float, default: null):
一个子采样因子，范围是（0，1）。如果指定的话，不能指定 "numPixels"。默认为不进行子采样。
numPixels（长，默认为空）。
要采样的像素的大概数量。如果指定，必须不指定'因子'。
seed（整数，默认为0）。
用于子采样的随机化种子。
dropNulls（布尔值，默认为true）。
对结果进行后置过滤，放弃具有空值属性的特征。
tileScale（浮点数，默认为1）。
用于减少聚合瓦片大小的比例因子；使用较大的瓦片尺度（例如2或4）可能会使计算在默认情况下耗尽内存。
geometries（布尔值，默认：false）。
如果为真，则添加采样像素的中心作为输出特征的几何属性。否则，geometries将被省略（节省内存）。
返回。特征集合
remap(from, to, defaultValue, bandName)
从输入值到输出值的映射，由两个平行列表表示。任何不包括在输入列表中的输入值，如果给出了defaultValue，则被设置为defaultValue，如果没有，则被屏蔽。请注意，由于浮点精度错误，包含浮点值的输入有时可能无法匹配。
参数。
this:image（图像）。
应用重映射的图像。
from (列表)。
源值（数字或ee.Array）。这个列表中的所有值将被映射到'to'中的相应值。
to (列表):
目标值（数字或ee.Array）。这些被用来替换'from'中的相应值。必须有与'from'相同数量的值。
defaultValue（对象，默认：null）。
用于替换'from'中没有匹配的值的默认值。如果不指定，未匹配的值会被屏蔽掉。
bandName（字符串，默认：null）。
要重新映射的乐队的名称。如果不指定，则使用图像中的第一个波段。
返回。图像
代码：
var geometry_roi = /* color: #d63000 *//* shown: false */ee.Geometry.Polygon([[[101.00841298710328, 3.5840956245106925],[101.15672841679078, 3.3538056172134443],[101.30779042850953, 3.4442833673276585],[101.11003652225953, 3.655364290850406],[101.15672841679078, 3.710182494184478],[101.02763906132203, 3.8061061368075664],[100.87932363163453, 3.6992191239822603]]]);// Define ROI
var roi=geometry_roi;
Map.centerObject(roi,10);//Date
var startDate = ee.Date('2019-01-01');
var endDate = ee.Date('2020-12-31');//S2影像的设定
// 为2019-2020年期间创建S2图像集。
var S2 = ee.ImageCollection('COPERNICUS/S2')//filter start and end date.filter(ee.Filter.date(startDate, endDate)).filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than',100)//filter according to drawn boundary.filterBounds(roi)print(S2.limit(10))
print(S2.aggregate_array('SPACECRAFT_NAME'))
// Function to calculate and add an NDVI band
var addNDVI = function(image) {return image.addBands(image.normalizedDifference(['B8', 'B4'] )); //'B8', 'B4'
}; // Add NDVI band to image collection
var S2 = S2.map(addNDVI).select(['nd']);
print('S2',S2.limit(10)) ;
var NDVI=S2.select('nd');// For month
var month = 1;// 计算区间数，这里用了一个很简单的方法进行了时间的获取
var months = endDate.difference(startDate,'month').divide(month).toInt();
//建立时间序列
var sequence = ee.List.sequence(0, months);
print(sequence)var sequence_s1 = sequence.map(function(num){
//数字化其变量num = ee.Number(num);
//中间间隔的时间设定，开始时间和就是默认的时间，终止时间往后加一个月var Start_interval = startDate.advance(num.multiply(month), 'month');var End_interval = startDate.advance(num.add(1).multiply(month), 'month');
//按照上面的时间进行影像筛选var subset = NDVI.filterDate(Start_interval,End_interval);
//获取子集的最大值，然后设定时间属性，因为在计算了max之后时间贤惠能干会丢失return subset.max().set('system:time_start',Start_interval);
});print('sequence_s1',sequence_s1)
//将上面的影像子集放入到影像集合中
var byMonthYear = ee.ImageCollection.fromImages(sequence_s1);print('byMonthYear',byMonthYear)
//这里有一个非常实用的函数将一个影像集合转化为一个单波段的影像
var multibandNDVI = byMonthYear.toBands().clip(roi);
print('multiband',multibandNDVI);//定义时间波段，将上面的子集命名的波段改名
var bandsName=['2019-01','2019-02','2019-03','2019-04','2019-05','2019-06','2019-07','2019-08','2019-09','2019-10','2019-11','2019-12','2020-01','2020-02','2020-03','2020-04','2020-05','2020-06','2020-07','2020-08','2020-09','2020-10','2020-11','2020-12']var multiband1_ndvi = multibandNDVI.rename(bandsName).clip(roi);//(monList)////s1影像的设定，基本步骤同上
var sentinel1_vh = ee.ImageCollection('COPERNICUS/S1_GRD').filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH')).select('VH').filter(ee.Filter.eq('instrumentMode', 'IW')).filter(ee.Filter.eq('resolution_meters', 10)).filter(ee.Filter.date(startDate, endDate)).filter(ee.Filter.bounds(roi))print('s1',sentinel1_vh);// For month
var month = 1;// Calculating number of intervals
var months = endDate.difference(startDate,'month').divide(month).toInt();
// Generating a sequence
var sequence = ee.List.sequence(0, months);
print(sequence)var sequence_s1 = sequence.map(function(num){num = ee.Number(num);var Start_interval = startDate.advance(num.multiply(month), 'month');var End_interval = startDate.advance(num.add(1).multiply(month), 'month');var subset = sentinel1_vh.filterDate(Start_interval,End_interval);return subset.median().set('system:time_start',Start_interval);
});print('sequence_s1',sequence_s1)
var byMonthYearS1 = ee.ImageCollection.fromImages(sequence_s1);
var multibands1 = byMonthYearS1.toBands().clip(roi);var multibands1 = multibands1.rename(bandsName).clip(roi);//.rename(monLists1).clip(roi);////将S1和S2合并
var combinedband=multiband1_ndvi.addBands(multibands1);
print('combinedband',combinedband);//将影像进行采样var training = combinedband.sample({region: roi,scale: 10,numPixels: 3000,tileScale:8,// geometries:true});//Map.addLayer(training,{},'points')//这里用Kmeans分类器进行训练
var clusterer = ee.Clusterer.wekaKMeans(20).train({features:training});//使用经过训练的聚类器对输入进行聚类。
var result_cluster =combinedband.cluster(clusterer).byte();//combands
// print('result_s2',result_s2)var clusters = [0, 1, 2, 3, 4,5,6,7,8,9,
10,11,12,13,14,15,16,17,18,19,
20,21,22,23,24,25,26,27,28,29];var values0 = [1, 2, 3, 4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30];var values1 = [0, 0, 0, 1,0,1,0,1,0,0,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0];
/*
var values1 = [0, 0, 0, 0,1,0,1,0,1,0,0,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0]; */
var values2 = [0, 0, 0, 2,0,3,0,4,0,0,0,0,1,0,0,2,3,0,0,0,0,0,0,4,0,0,0,0,0,0];
/*
var values2 = [0, 0, 0, 0, 2,0,3,0,4,0,0,0,0,1,0,0,2,3,0,0,0,0,0,0,4,0,0,0,0,0];*/ //这里根据定义的value值进行重新赋值映射，
//var result_cluster = result_cluster.remap(clusters,values0);
var remapped_cluster = result_cluster.rename('remapped');// .remap(clusters, values0).byte().clip(roi);//.updateMask(slope.lt(slope_th))
var remapped_cluster1=remapped_cluster.remap(values0,values1);
var remapped_cluster1=remapped_cluster1.updateMask(remapped_cluster1);
var remapped_cluster2=remapped_cluster.remap(values0,values2).updateMask(remapped_cluster1);
/////
var comb_ndvi_cluster=multiband1_ndvi.addBands(remapped_cluster);
print("comb_ndvi_cluster",comb_ndvi_cluster);
// 定义图表的定制选项。
var options = {lineWidth: 1,pointSize: 2,hAxis: {title: 'Year-Month'},vAxis: {title: 'NDVI'},title: 'Sentinel-2 NDVI spectra in classified regions'
};// 制作图表，设置选项。
var chart_class_ndvi = ui.Chart.image.byClass(comb_ndvi_cluster, 'remapped', roi, ee.Reducer.median(), 500)//, classNames, wavelengths).setOptions(options).setChartType('ScatterChart');// Print the chart.
//print(chart_class_ndvi);//
var comb_vh_cluster=multibands1.addBands(remapped_cluster);
print("comb_vh_cluster",comb_vh_cluster);
// Define chart customization options.
var options = {lineWidth: 1,pointSize: 2,hAxis: {title: 'Year-Month'},vAxis: {title: 'VH'},title: 'Senetinel-1 VH spectra in classified regions'
};// Make the chart, set the options.
var chart_class_vh = ui.Chart.image.byClass(comb_vh_cluster, 'remapped', roi, ee.Reducer.median(), 500)//, classNames, wavelengths).setOptions(options).setChartType('ScatterChart');// Print the chart.
//print(chart_class_vh);Map.addLayer(multiband1_ndvi , {min: 0.2, max: 0.8}, 'NDVI',0);
Map.addLayer(multibands1 , {min: -25, max: -10}, 'VH',0);Map.addLayer(remapped_cluster.randomVisualizer(), {}, 're_groups_s2',0);
Map.addLayer(remapped_cluster1, {min:1,max:1,palette:['white','red']}, 'binary_rice',1);
Map.addLayer(remapped_cluster2.randomVisualizer(), {}, 'pattern_rice',1); ///// 在地图旁边创建一个面板
var panel = ui.Panel({style:{width: '550px'}});
ui.root.add(panel);// value: 'SENTINEL-2 NDVI & SENTINEL-1 VH Explorer',var label1 = ui.Label({value: 'SENTINEL-2 NDVI & SENTINEL-1 VH Explorer: Paddy Clustering',style: {//fontSize: '20px',color: 'white',backgroundColor :'blue',// border: '1px solid black',// fontWeight: 'bold',padding: '5px',margin: '12px 0px 0px 8px', fontWeight: 'bold', fontSize: '14px', border: '1px solid black', padding: '5px 5px 5px 5px',}
// style: {// fontSize: '16px',// color: 'red',// fontWeight: 'bold',// padding: '5px',// }
});var label2 = ui.Label({value: 'Developed by Rudiyanto',style: {fontSize: '16px',color: 'red',// fontWeight: 'bold',padding: '5px',}
});var affiliationLabel= ui.Label('Instructions',{fontWeight: 'bold',padding: '3px',});
var affiliation= ui.Label("Program of Crop Science\n"+"Faculty of Fisheries and Food Science\n"+"Universiti Malaysia Terengganu\n"+"Email: rudiyanto@umt.edu.my\n", {whiteSpace:'pre',padding: '5px',fontSize: '16px',color: 'red',// fontWeight: 'bold',});var affiliationPanel = ui.Panel([affiliationLabel,affiliation]);
var affiliationPanel = ui.Panel([affiliation]);panel.add(label1);panel.add(label2);panel.add(affiliationPanel);panel.add(chart_class_ndvi);panel.add(chart_class_vh);这是将两个影像合并在一起并改名的波段：
S1影像
S2影像
整体的代码界面：
结论
- 由于使用了 Sentinel-1 和 2 时间序列数据，它不仅可以生成高分辨率的水稻分布图，还可以绘制热带地区的强度和种植日历；
- 由于在 GEE 中使用了无监督的 K-Means 方法，它可以在没有大量实地调查数据的情况下快速生成高精度地图；
- 由于使用了感兴趣的目标区域，它减少了盐和胡椒的影响；
- 由于使用了每月最绿的图像或 Sentinel-2 数据的最大 NDVI 值组合，因此无需遮盖云层和云层阴影；
- 由于 Sentinel-1 和 2 时间序列数据的融合，它能够识别水稻物候阶段的特征，包括土壤耕作和种植、营养、生殖和成熟阶段。
文献引用：
Fatchurrachman; Rudiyanto; Soh, N.C.; Shah, R.M.; Giap, S.G.E.; Setiawan, B.I.; Minasny, B. High-Resolution Mapping of Paddy Rice Extent and Growth Stages across Peninsular Malaysia Using a Fusion of Sentinel-1 and 2 Time Series Data in Google Earth Engine. Remote Sens. 2022, 14, 1875. https://doi.org/10.3390/rs14081875
Google Earth Engine（GEE）——Sentinel-1 和 2 数据的融合，水稻范围识别和水稻种植季节区分地图绘制—马来西亚为例相关推荐
- 使用GEE（Google Earth Engine）下载STRM30m分辨率高程数据
简介 本文介绍如何使用GEE(Google Earth Engine)下载STRM30m分辨率高程数据并用ArcMap加载行政边界裁剪栅格图层. 准备工作 申请Google账号 注册Google Dr ...
- Google Earth Engine(GEE)——User memory limit exceeded（2）
上一次我们已经知道如何去进行避免这种错误的发生,有关详细内容,如果单单只是解决这个问题我们用到的是limit 和 first,上一次的博客在这里: (207条消息) Google Earth Engi ...
- Google Earth Engine(GEE) 01-中输入提示快捷键Ctrl+space无法使用的问题
Google Earth Engine(GEE) 01-中输入提示快捷键Ctrl+space无法使用的问题 GEE中 Ctrl+space组合键用于代码输入快捷提示,能够提高编码的准确度和速度,但是, ...
- Google Earth Engine(GEE)批量下载代码（以 NDVI数据为例）
下载数据先准备工作(具体细节都能查到):科学上网 谷歌邮箱,谷歌邮箱注册GEE账号. 一:导入需要下载边界shp文件. 标题 找到自己的shp文件,导入除了sbx文件的所有文件. 导入成功 命名ass ...
- 使用Google Earth Engine (GEE)实现MODIS数据批量下载
使用Google Earth Engine GEE实现MODIS数据批量下载 前言 下载数据代码 批量执行run任务 关注公众号,分享GIS知识.ArcGIS教程.SCI论文与科研日常等 前言 上图是 ...
- 基于google earth engine(GEE)下载研究区域影像
基于google earth engine(GEE)下载研究区域影像 当研究需要Landsat数据时,我们可以通过USGS官网或者地理空间数据云平台下载.由于地理空间数据云目前无法下载到较新的数据,可 ...
- Google Earth Engine(GEE)——可视化动态图
代码: var geometry = /* color: #d63000 *//* shown: false *//* displayProperties: [{"type": & ...
- Google Earth Engine (GEE) ——卫星影像的监督分类（svm）
问题 GEE 提供哪些机器学习技术? 如何对卫星图像进行监督分类? 如何评估分类器的准确性? 如何手动创建自己的几何图形? 目标 练习查找无云图像和使用手绘几何导入 学习训练和应用分类算法所需的基本功 ...
- 关于google earth engine(GEE)的一些想法与大胆预测
我接触GEE有两年了,GEE留给我的印象是:无所不能. 不管是从庞大的数据量,还是包含遥感的各类算法:随机森林.SVM.CNN,都让人惊讶. 从GEE的云端操作来看,传统遥感需要几个月做出来的全国ND ...
最新文章
- JAVA偏向锁的什么时候释放_Java中的偏向锁
- 一致性哈希算法——算法解决的核心问题是当slot数发生变化时，能够尽量少的移动数据
- fastjson 添加key value_Redis中设置了过期时间的Key，那么你还要知道些什么？
- C++容器（三）：pair类型
- hive遇到FAILED: Execution Error, return code 2 from org.apache.hadoop.hive.ql.exec.mr.MapRedTask错误...
- java redis 原子操作_7000字 Redis 超详细总结笔记 ！建议收藏
- vba判断是否为数字的方法小集
- C++ vector和list的区别
- 自动编译 打包 签名 android程序
- SQL执行insert into后返回主键ID
- C#表驱动法+一点反射实现“得到指定位数随机不重复字符串”三种方式的封装...
- 微波遥感SNAP（一）——基于Sentinel-1雷达数据重建数字高程模型
- Spring框架RCE 0day漏洞，官方解决方案
- css宋体代码_css怎么设置字体为宋体
- 实验一 路由器的基本管理
- 易语言-通过快捷方式.lnk文件获取.exe源文件地址
- 17个国内外手机APP应用开发平台推荐
- python3.5项目升级3.8问题汇总
- 无约束低分辨率人脸识别综述一：用于低分辨率人脸识别的数据集
- 主机和服务器之间有何区别？