Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析

2024-08-31 10:36

本文主要是介绍Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

简介

函数

ee.Array.identity(size)

Arguments:

Returns: Array

transpose(axis1, axis2)

Arguments:

Returns: Array

matrixMultiply(image2)

Arguments:

Returns: Image

matrixSolve(image2)

Arguments:

Returns: Image

arrayFlatten(coordinateLabels, separator)

Arguments:

Returns: Image

arrayReduce(reducer, axes, fieldAxis)

Arguments:

Returns: Image

代码

结果


简介

惠特克(GEE)平滑算法是一种用于时间序列预测的统计方法,特别适用于非线性、非平稳和非高斯的数据。该算法基于广义估计方程,通过最小化残差的平方和来拟合数据并找到最佳的平滑曲线。

GEE平滑算法的主要思想是在时间序列数据中引入一个平滑函数来描述数据的趋势和周期性变化。该平滑函数由一系列基函数的线性组合组成,其中每个基函数具有不同的频率和振幅。通过调整基函数的权重,可以得到最佳的平滑曲线,以最大程度地拟合数据。

在实际应用中,GEE平滑算法通常与其他统计方法结合使用,例如自回归移动平均模型(ARIMA)或指数平滑法。通过将GEE平滑算法与其他方法相结合,可以进一步提高时间序列的预测准确度和稳定性。

总的来说,GEE平滑算法是一种针对非线性、非平稳和非高斯数据的时间序列预测方法,通过引入一个平滑函数来描述数据的趋势和周期性变化,以最大程度地拟合数据。它在实际应用中通常与其他统计方法结合使用,以进一步提高预测的准确度和稳定性。

函数

ee.Array.identity(size)

Creates a 2D identity matrix of the given size.

创建一个给定大小的二维标识矩阵。

Arguments:

size (Integer):

The length of each axis.

Returns: Array

transpose(axis1axis2)

Transposes two dimensions of an array.

平移数组的两个维度。

Arguments:

this:array (Array):

Array to transpose.

axis1 (Integer, default: 0):

First axis to swap.

axis2 (Integer, default: 1):

Second axis to swap.

Returns: Array

matrixMultiply(image2)

Returns the matrix multiplication A * B for each matched pair of bands in image1 and image2. If either image1 or image2 has only 1 band, then it is used against all the bands in the other image. If the images have the same number of bands, but not the same names, they're used pairwise in the natural order. The output bands are named for the longer of the two inputs, or if they're equal in length, in image1's order. The type of the output pixels is the union of the input types.

返回图像 1 和图像 2 中每对匹配波段的矩阵乘法 A * B。如果图像 1 或图像 2 中只有一个波段,则该波段将与另一幅图像中的所有波段相对应。如果图像中的条带数量相同,但名称不同,则按自然顺序成对使用。输出波段以两个输入波段中较长的一个命名,如果两个输入波段长度相等,则按图像 1 的顺序命名。输出像素的类型是输入类型的组合。

Arguments:

this:image1 (Image):

The image from which the left operand bands are taken.

image2 (Image):

The image from which the right operand bands are taken.

Returns: Image

matrixSolve(image2)

Solves for x in the matrix equation A * x = B, finding a least-squares solution if A is overdetermined for each matched pair of bands in image1 and image2. If either image1 or image2 has only 1 band, then it is used against all the bands in the other image. If the images have the same number of bands, but not the same names, they're used pairwise in the natural order. The output bands are named for the longer of the two inputs, or if they're equal in length, in image1's order. The type of the output pixels is the union of the input types.

求解矩阵方程 A * x = B 中的 x,如果 A 对图像 1 和图像 2 中每对匹配的波段都是过确定的,则找到最小二乘法解。如果图像 1 或图像 2 中只有一个波段,则使用该波段与另一幅图像中的所有波段进行比对。如果图像中的波段数相同,但名称不相同,则按自然顺序成对使用。输出波段以两个输入波段中较长的一个命名,如果两个输入波段长度相等,则按图像 1 的顺序命名。输出像素的类型是输入类型的组合。

Arguments:

this:image1 (Image):

The image from which the left operand bands are taken.

image2 (Image):

The image from which the right operand bands are taken.

Returns: Image

arrayFlatten(coordinateLabels, separator)

Converts a single-band image of equal-shape multidimensional pixels to an image of scalar pixels, with one band for each element of the array.

将等形多维像素的单波段图像转换为标量像素图像,阵列中的每个元素对应一个波段。

Arguments:

this:image (Image):

Image of multidimensional pixels to flatten.

coordinateLabels (List):

Name of each position along each axis. For example, 2x2 arrays with axes meaning 'day' and 'color' could have labels like [['monday', 'tuesday'], ['red', 'green']], resulting in band names'monday_red', 'monday_green', 'tuesday_red', and 'tuesday_green'.

separator (String, default: "_"):

Separator between array labels in each band name.

Returns: Image

arrayReduce(reducer, axes, fieldAxis)

Reduces elements of each array pixel.

减少每个阵列像素的元素。

Arguments:

this:input (Image):

Input image.

reducer (Reducer):

The reducer to apply.

axes (List):

The list of array axes to reduce in each pixel. The output will have a length of 1 in all these axes.

fieldAxis (Integer, default: null):

The axis for the reducer's input and output fields. Only required if the reducer has multiple inputs or outputs.

Returns: Image

代码

//加载研究区
var geometry = /* color: #d63000 *//* displayProperties: [{"type": "rectangle"}] */ee.Geometry.Polygon([[[113.44773227683973, 38.6708907304602],[113.44773227683973, 38.64783484482313],[113.47588474266004, 38.64783484482313],[113.47588474266004, 38.6708907304602]]], null, false);// 将 qa 位图像转换为标志的辅助函数
function extractBits(image, start, end, newName) {// 计算我们需要提取的比特。var pattern = 0;for (var i = start; i <= end; i++) {pattern += Math.pow(2, i);}// 返回提取的质量保证位的单波段图像,并为该波段命名。return image.select([0], [newName]).bitwiseAnd(pattern).rightShift(start);
}// 在输入矩阵上获取指定阶次的差分矩阵的函数。将矩阵和阶次作为参数
function getDifferenceMatrix(inputMatrix, order){var rowCount = ee.Number(inputMatrix.length().get([0]));var left = inputMatrix.slice(0,0,rowCount.subtract(1));var right = inputMatrix.slice(0,1,rowCount);if (order > 1 ){return getDifferenceMatrix(left.subtract(right), order-1)}return left.subtract(right);
};// 将数组图像解包为图像和波段
// 以数组图像、图像 ID 列表和乐队名称列表为参数
function unpack(arrayImage, imageIds, bands){function iter(item, icoll){function innerIter(innerItem, innerList){return ee.List(innerList).add(ee.String(item).cat("_").cat(ee.String(innerItem)))}var temp = bands.iterate(innerIter, ee.List([]));return ee.ImageCollection(icoll).merge(ee.ImageCollection(ee.Image(arrayImage).select(temp,bands).set("id",item)))}var imgcoll  = ee.ImageCollection(imageIds.iterate(iter, ee.ImageCollection([])));return imgcoll}// 用于计算回归结果的反对数比率并转换回百分比单位的函数
function inverseLogRatio(image) {var bands = image.bandNames();var t = image.get("system:time_start");var ilrImage = ee.Image(100).divide(ee.Image(1).add(image.exp())).rename(bands);return ilrImage.set("system:time_start",t);
}function whittakerSmoothing(imageCollection, isCompositional, lambda){// 快速配置以设置默认值if (isCompositional === undefined || isCompositional !==true) isCompositional = false;if (lambda === undefined ) lambda = 10;// 程序启动  var ic = imageCollection.map(function(image){var t = image.get("system:time_start");return image.toFloat().set("system:time_start",t);});var dimension = ic.size();var identity_mat = ee.Array.identity(dimension);var difference_mat = getDifferenceMatrix(identity_mat,3);var difference_mat_transpose = difference_mat.transpose();var lamda_difference_mat = difference_mat_transpose.multiply(lambda);var res_mat = lamda_difference_mat.matrixMultiply(difference_mat);var hat_matrix = res_mat.add(identity_mat);// 备份原始数据var original = ic;// 获取原始图像属性var properties = ee.List(ic.iterate(function(image, list){return ee.List(list).add(image.toDictionary());},[]));var time = ee.List(ic.iterate(function(image, list){return ee.List(list).add(image.get("system:time_start"));},[]));// 如果数据是合成的// 计算图像在 0 到 100 之间的对比率。首先// 夹在 delta 和 100-delta 之间,其中 delta 是一个很小的正值。if (isCompositional){ic = ic.map(function(image){var t = image.get("system:time_start");var delta = 0.001;var bands = image.bandNames();image = image.clamp(delta,100-delta);image = (ee.Image.constant(100).subtract(image)).divide(image).log().rename(bands);return image.set("system:time_start",t);});}var arrayImage = original.toArray();var coeffimage = ee.Image(hat_matrix);var smoothImage = coeffimage.matrixSolve(arrayImage);var idlist = ee.List(ic.iterate(function(image, list){return ee.List(list).add(image.id());},[]));var bandlist = ee.Image(ic.first()).bandNames();var flatImage = smoothImage.arrayFlatten([idlist,bandlist]);var smoothCollection = ee.ImageCollection(unpack(flatImage, idlist, bandlist));if (isCompositional){smoothCollection = smoothCollection.map(inverseLogRatio);}// 通过添加后缀fitted获得新的乐队名称var newBandNames = bandlist.map(function(band){return ee.String(band).cat("_fitted")});// 重新命名平滑图像中的波段smoothCollection = smoothCollection.map(function(image){return ee.Image(image).rename(newBandNames)});// 一个非常笨的方法,可以flatten谷歌地球引擎生成的 ID,这样就可以将两张图片合并为图表了var dumbimg = arrayImage.arrayFlatten([idlist,bandlist]);var dumbcoll = ee.ImageCollection(unpack(dumbimg,idlist, bandlist));var outCollection = dumbcoll.combine(smoothCollection);var outCollectionProp = outCollection.iterate(function(image,list){var t = image.get("system:time_start")return ee.List(list).add(image.set(properties.get(ee.List(list).size())));},[]);var outCollectionProp = outCollection.iterate(function(image,list){return ee.List(list).add(image.set("system:time_start",time.get(ee.List(list).size())));},[]);var residue_sq = smoothImage.subtract(arrayImage).pow(ee.Image(2)).divide(dimension);var rmse_array = residue_sq.arrayReduce(ee.Reducer.sum(),[0]).pow(ee.Image(1/2));var rmseImage = rmse_array.arrayFlatten([["rmse"],bandlist]);return [ee.ImageCollection.fromImages(outCollectionProp), rmseImage];
}var ndvi =ee.ImageCollection("NOAA/VIIRS/001/VNP13A1").select('NDVI').filterDate("2019-01-01","2019-12-31");
// 去除屏蔽像素
ndvi = ndvi.map(function(img){return img.unmask(ndvi.mean())});var ndvi =  whittakerSmoothing(ndvi)[0];// 添加图表
print(ui.Chart.image.series(ndvi.select(['NDVI', 'NDVI_fitted']), geometry, ee.Reducer.mean(), 500).setSeriesNames(['NDVI', 'NDVI_fitted']).setOptions({title: 'smoothed',lineWidth: 1,pointSize: 3,
}));

结果

这篇关于Google Earth Engine:对NDVI进行惠特克平滑算法进行长时序分析的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1123663

相关文章

Nginx中配置使用非默认80端口进行服务的完整指南

《Nginx中配置使用非默认80端口进行服务的完整指南》在实际生产环境中,我们经常需要将Nginx配置在其他端口上运行,本文将详细介绍如何在Nginx中配置使用非默认端口进行服务,希望对大家有所帮助... 目录一、为什么需要使用非默认端口二、配置Nginx使用非默认端口的基本方法2.1 修改listen指令

MySQL按时间维度对亿级数据表进行平滑分表

《MySQL按时间维度对亿级数据表进行平滑分表》本文将以一个真实的4亿数据表分表案例为基础,详细介绍如何在不影响线上业务的情况下,完成按时间维度分表的完整过程,感兴趣的小伙伴可以了解一下... 目录引言一、为什么我们需要分表1.1 单表数据量过大的问题1.2 分表方案选型二、分表前的准备工作2.1 数据评估

MySQL进行分片合并的实现步骤

《MySQL进行分片合并的实现步骤》分片合并是指在分布式数据库系统中,将不同分片上的查询结果进行整合,以获得完整的查询结果,下面就来具体介绍一下,感兴趣的可以了解一下... 目录环境准备项目依赖数据源配置分片上下文分片查询和合并代码实现1. 查询单条记录2. 跨分片查询和合并测试结论分片合并(Shardin

Android 缓存日志Logcat导出与分析最佳实践

《Android缓存日志Logcat导出与分析最佳实践》本文全面介绍AndroidLogcat缓存日志的导出与分析方法,涵盖按进程、缓冲区类型及日志级别过滤,自动化工具使用,常见问题解决方案和最佳实... 目录android 缓存日志(Logcat)导出与分析全攻略为什么要导出缓存日志?按需过滤导出1. 按

Linux中的HTTPS协议原理分析

《Linux中的HTTPS协议原理分析》文章解释了HTTPS的必要性:HTTP明文传输易被篡改和劫持,HTTPS通过非对称加密协商对称密钥、CA证书认证和混合加密机制,有效防范中间人攻击,保障通信安全... 目录一、什么是加密和解密?二、为什么需要加密?三、常见的加密方式3.1 对称加密3.2非对称加密四、

MySQL中读写分离方案对比分析与选型建议

《MySQL中读写分离方案对比分析与选型建议》MySQL读写分离是提升数据库可用性和性能的常见手段,本文将围绕现实生产环境中常见的几种读写分离模式进行系统对比,希望对大家有所帮助... 目录一、问题背景介绍二、多种解决方案对比2.1 原生mysql主从复制2.2 Proxy层中间件:ProxySQL2.3

SpringBoot结合Knife4j进行API分组授权管理配置详解

《SpringBoot结合Knife4j进行API分组授权管理配置详解》在现代的微服务架构中,API文档和授权管理是不可或缺的一部分,本文将介绍如何在SpringBoot应用中集成Knife4j,并进... 目录环境准备配置 Swagger配置 Swagger OpenAPI自定义 Swagger UI 底

基于Python Playwright进行前端性能测试的脚本实现

《基于PythonPlaywright进行前端性能测试的脚本实现》在当今Web应用开发中,性能优化是提升用户体验的关键因素之一,本文将介绍如何使用Playwright构建一个自动化性能测试工具,希望... 目录引言工具概述整体架构核心实现解析1. 浏览器初始化2. 性能数据收集3. 资源分析4. 关键性能指

Nginx进行平滑升级的实战指南(不中断服务版本更新)

《Nginx进行平滑升级的实战指南(不中断服务版本更新)》Nginx的平滑升级(也称为热升级)是一种在不停止服务的情况下更新Nginx版本或添加模块的方法,这种升级方式确保了服务的高可用性,避免了因升... 目录一.下载并编译新版Nginx1.下载解压2.编译二.替换可执行文件,并平滑升级1.替换可执行文件

python使用Akshare与Streamlit实现股票估值分析教程(图文代码)

《python使用Akshare与Streamlit实现股票估值分析教程(图文代码)》入职测试中的一道题,要求:从Akshare下载某一个股票近十年的财务报表包括,资产负债表,利润表,现金流量表,保存... 目录一、前言二、核心知识点梳理1、Akshare数据获取2、Pandas数据处理3、Matplotl