soderayer的博客

Python遥感图像处理主要采用python编程语言，使用arcpy以及gdal等库进行遥感图像数据处理及应用。包括遥感图像的拼接、镶嵌、裁剪、格式转换等基础与处理操作以及常用遥感指数的批量计算，比如NDVI\EVI等遥感指数的计算与应用。遥感影像数据往往数量众多，而且是具有标准分幅大小的，具有统一的格式、坐标系统等标准化信息，因此，通常需要对遥感图像进行大批量处理。本专栏主要采用批处理的方式去完成。比如使用数据为MRT处理后的tif格式Modis数据，前期批下载、批处理具体处理过程之前博文有所介绍

完整的遥感图像滤波分析通常包括以下几个步骤和方法：噪声检测：首先需要检测图像中存在的噪声类型，例如高斯噪声、椒盐噪声等。可以使用统计分析或图像处理技术进行噪声检测。噪声去除：根据噪声类型选择合适的去噪方法，例如：高斯滤波器（Gaussian Filter）：适用于高斯噪声，通过应用高斯函数对图像进行平滑。中值滤波器（Median Filter）：适用于椒盐噪声，通过取邻域内像素值的中值来平滑图像。自适应滤波器（Adaptive Filter）：通过动态调整滤波参数，根据图像的局部特性进行去噪。

遥感图像的统计特征是对图像中像素值的统计分布进行定量化描述的过程。这些统计特征可以提供关于图像内容和特性的有用信息。下面是一些常用的遥感图像统计特征描述方法：平均值（Mean）：计算图像中所有像素值的平均值，可以反映整个图像的亮度水平。方差（Variance）：计算图像中所有像素值与其平均值之间的偏差的平方的平均值，用于描述图像的对比度和纹理变化程度。标准差（Standard Deviation）：方差的平方根，衡量像素值在平均值附近的分散程度。

使用多波段遥感图像进行主成分分析，这里使用了6个波段的数据计算和显示效果图如下：左边是原图像IR+R+G显示图，右边是计算得到的3个主成分组合显示的RGB图像。

方向梯度直方图HOG（Histogram of Oriented Gradient, HOG）特征是一种在计算机视觉和图像处理中用来进行物体检测的特征描述子。它通过计算和统计图像局部区域的梯度方向直方图来构成特征。Hog特征结合SVM分类器已经被广泛应用于图像识别中，尤其在行人检测中获得了极大的成功。与其他的特征描述方法相比，HOG有很多优点。首先，由于HOG是在图像的局部方格单元上操作，所以它对图像几何的和光学的形变都能保持很好的不变性，这两种形变只会出现在更大的空间领域上。

局部二值模式（Local Binary Pattern, LBP）是一种用于描述图像纹理特征的算法。它通过对图像的每个像素点与其邻域像素进行比较，得到一个二进制编码来表示该像素点的纹理信息。LBP 算法的优点是简单高效，并且对噪声有一定的鲁棒性。它广泛应用于图像纹理分析、人脸识别、行人检测、纹理分类等领域。使用 LBP 特征可以准确地描述图像纹理信息，提取出图像中重要的纹理特征，为后续的图像处理和分析任务提供有用的信息。

在遥感图像处理中，计算纹理特征是一种常见的方法来描述图像中不同区域之间的纹理差异。本文主要是用GLCM来计算遥感图像的纹理特征。

部分原文介绍如下，我们可以通过一下参数来计算地表反射率以及进行地表温度单位转换。温度产品展示图（上文中链接可查看）Scaling Factor是低保温度计算的关键参数，上表中分别介绍了Collection1和Collection 2两种数据的缩放比例，两种数据使用了不同的缩放因子，使用的是开尔文温度。开尔文（Kelvins），为热力学温标或称绝对温标，是中的温度单位。开尔文温度和华氏温度转换公式如下：从开氏温标换算至其他温度单位从其他温度单位换算至开氏温标摄氏温标华氏温标兰金温标。

GF-2数据，该数据经过辐射定标、大气校正、影像融合和正射校正处理。

运行环境运行环境：windows10 pycharm python3.7.7包等。实验目标实现遥感图像实现邻域分析中的最多数值分析，设计背景值区域的像元不参与计算。

根据该函数的参数，要实现图像的平移，其实我们只需要修改图像的左上角起始点坐标就能够实现了，比如我们将图像向东易懂5000m，向南移动5000m，那么，我们是需要将图像六元组中的第一个和第四个要素平移一定距离就行了。Y坐标方向是负值，所以我们向南移动就是-5000m，向东移动就是+5000m。元组是不可变的，所以我们需要将元组转换为列表，更改列表，然后将列表转换回元组。运行环境：windows10 pycharm python3.7.7。GDAL 中，GeoTransform是一个六个元素的元组。

最近在调试一段代码，该代码输入数据为遥感图像，但是计算出的图像分类结果使用jpg（png等）图片格式保存的输出结果，作为一个遥感工作者，显然更想得到一个具有空间位置信息的遥感图像，比如Geotiff格式。另一种方法是直接将Jpg格式的图片转化为带坐标的遥感图像。这里需要使用到一个参考影像数据（可以将原始输入计算的遥感数据作为参考数据，主要用于获取待保存数据的空间信息。1.Jpg格式结果具有多种类别结果，每种类别具有不同的颜色，所以首先我们将图像灰度化，可以使用opencv库来实现，将转化后的图像保存下来。

将原始遥感图像各个波段拉伸到0-255范围显示。运行环境：windows10 pycharm python3.7.7包等。

计算思路：（1）分别提取参考影像和分类影像像元值，通过GDAL读取栅格影像数据的像元值，使用pandas保存到dataframe中，第一列使用FID记录像元ID,第二列记录像元值。（2）合并两个dataframe为新的dataframe，参考影像中像元值包含实际类别以及背景值0值，我们需要删除背景值所对应的像元，因此，根据删除像元值为0的数据。（3）分别读取上述结果中的参考像元和分类像元对应的值，获得参考值和分类值两个序列的数值。（4）对上述两个系列的值进行精度计算和制图可视化分析。

续-上一篇中使用csv文件计算的相关系数热力图，本篇我们直接使用遥感图像来计算图像波段之间的相关系数。

批量读取原始压缩包解压后的数据，然后进行波段组合运算。文件列表如下：之前写过的是基于Landsat8数据，将所有tif文件放置在一个文件夹种实现的波段组合。文件路径很多，我们需要先读取到需要组合波段的完整路径。选取以上9个波段进行波段组合。最后读取到的是全部的波段完整路径，部分显示如下：接着就可以开始按照上一篇文章之中的代码进行波段组合运算了。

本篇对上一篇博客代码进行少量修改，实现遥感图像按照指定大小等距裁剪。比如我们需要将原始遥感图像按照60km*60km大小标准块裁剪，我们将上一篇中的裁剪函数参数修改一下，将行列号修改为裁剪大小。

之前使用arcpy python2.7写了一篇进行遥感数据波段组合计算的博文，使用场景是将所有单波段数据放在一个文件夹中（如下图），文件名称前缀一样，后缀波段编号用来区分不同的波段，数据按照名称默认排序，读取的时候也是按照排列顺序读取。

给多光谱遥感图像各个波段绘制热力图，首先需要计算波段之间的相关系数矩阵，而计算遥感图像波段相关系数矩阵有不同的方法，常用的我们可以采用遥感图像处理软件计算，比如ENVI软件就可以计算相关系数矩阵，使用工具箱中的Statistics工具即可进行多种统计运算。我们这里不使用遥感图像软件计算，直接使用Python计算。要实现整个目标，我们首先需要提取遥感图像波段中各个像元的信息。这里我们可以将遥感图像读取并保存为csv文件，实现这一步我们可以参照之前的一篇文章，链接如下：Python遥感图像处理应用篇(十四)：G

遥感影像拼接分多种不同的情况，比如比较常用的是两幅影像拼接或者多个图像合并为一个图像。GDAL中实现影像拼接的方式也有多种，比如常用可以通过numpy读取影像数组，并计算各个图像的numpy数组范围，合并多个numpy数组为一个最终的结果，并将其重新输出为一幅新的图像，即可实现图像的拼接。另外比较方便的一种方法便是使用GDAL.warp()函数方法。这个函数的功能很多，可以实现多种不同的目的，比如用于坐标系转换、投影变换、图像合并与镶嵌、地理范围裁剪、更改分辨率、矢量裁剪等多种功能。这些功能可以通过设置Op

给定一幅图像，我们需要将图像等距离分割裁剪，分割行列号相同，比如分割裁剪为N行N列的多幅遥感图像。例如原图如下： 裁剪为三行三列之后如下：设置不同波段组合为不同颜色以便区分边界。3*3裁剪图像运行环境：运行环境：运行环境：windows10 pycharm python3.7.7 GDAL-3.2.3-cp37-cp37m-win_amd64.whl包等。我们主要是用warp是扭曲或者变形的意思，该函数是对一个数据或者多个数据集进行变形，用法如下：参数包括：关键字参数是：如果选项是作为gdal.War

之前也写过NDVI批量计算的实现方法，采用的是Acrpy python2.7实现的，这里我们采用GDAL实现NDVI指数的批量计算。使用数据，多波段遥感图像数据。Landsat05 Colection2 Level2数据。NDVI=(NIR-RED)/(NIR+RED) 就算之前需要搞清楚多光谱数据中的波段对应关系。影像数据地表反射率修正系数：0.0000275+（-0.2）运行环境：运行环境：windows10 pycharm python3.7.7 GDAL-3.2.3-cp37-cp37m-w

1.实现目标批量将多波段图像提取为各个单波段图像运行环境：windows10 pycharmpython3.7.7GDAL-3.2.3-cp37-cp37m-win_amd64.whl包等。2.实现代码"""此段代码将遥感图像背景值去除处理"""import numpy as npfrom osgeo import gdal, gdalconstimport os#将遥感影像归一化处理 写成函数def GetEnvolopePoint(inputpath,outp...

最近下载了一些遥感影像数据，这些数据都包含大量的外围背景数据，如下图所示：外围背景值都为0值。本文描述的是采用python批量处理外围背景，只保留最小外围背景区域。如下图：基本思路如下：首先获取遥感影像有数据区域四个角的坐标，这里的坐标我们用行列号表示，即左上角（x1,y1),右上角（x2,y2），左下角（x3,y3），右下角（x4,y4）。我们需要通过遍历遥感数据的方法获取到四个角点的行列号。然后，根据四角坐标计算四至范围，即left，top，right，bottom值。之后，分别获取左上角坐标和右下角坐

所用输入数据为遥感影像数据，可以是单波段也可以是多波段数据。运行环境：windows10 pycharm python3.7.7 GDAL-3.2.3-cp37-cp37m-win_amd64.whl包等。完整代码如下：采用批处理方式，输入和输出为文件路径。只需要修改输入输出文件路径即可计算所需结果。3.数据运行测试与结果比如采用Landsat8中的NDVI数据进行批量计算. 图像对比：图像数值范围： 图像背景值变化：原图是nodata计算之后背景值出现变化，变为0值。 像元值查看对比：..

本文所使用数据和上一篇数据相同，上一篇链接：Python遥感图像处理应用篇(十六)：GDAL 将归一化处理csv数据转化为遥感影像_空中旋转篮球的博客-CSDN博客上一篇文章将csv文件转化为单波段图像，这一篇将csv文件合成为一个多波段图像。注意一个前提是需要知道原图像的基本信息，如数据行列号，投影信息，数据类型信息等。实现写入多光谱数据的时候，核心问题在于使用循环写入各个通道数据。这里直接写入的是二维数组数据，看别人的文章的时候，有的使用的三维数组，如格式（channel,row,col）,我在使用的时

将上一篇文章中得到的计算结果作为转换数据。链接如下：Python遥感图像处理应用篇(十五)：GDAL 读取多光谱数据做归一化处理_空中旋转篮球的博客-CSDN博客基本思路，先读取csv数据，现读取每一列数据为一维数组，然后根据原始影像数据的大小（row，col）将一维数组塑形，即使用reshape函数变换为（row，col）形状对的二维数组。使用gdal将二维数组转换为tif图像。完整代码如下：修改一下以上代码，可以将数据转化为和原始数据相同波段数的多波段遥感影像。可以使用for循环一个通道一个通道的

上一篇文中计算得到的csv文件：Python遥感图像处理应用篇(十四)：GDAL 读取多光谱数据为二维数组并存入csv文件_空中旋转篮球的博客-CSDN博客我们在做归一化处理的时候，遥感软件往往是针对单波段图像的，或者不方便批处理图像。我们可以使用python自己实现计算处理。 首先将数据读取为numpy数组，然后根据数组运算，计算得到归一化计算结果；这里采用（X-Xmin）/（Xmax-Xmin）计算公式归一化处理。计算完整代码如下： 2.3 批量计算代码将以上代码修改为批量计算代码，可以一次性

之前写过一遍文章是采用Arcgis读取多光谱数据为csv文件的，相对来说有一些麻烦，本文有个更简单一点的方法，使用python直接读取。本文使用Landsat05数据，采用GEE下载的多光谱数据，该数据下载下来包含19个波段。这些波段信息可以查看数据该数据详细说明。1-7是其对应的多光谱波段，第六波段是地表温度数据（ST),其他的是(SR) 地表反射率产品数据。.....................

1.绿度归一化植被指数GNDVI公式：GNDVI=TM4-TM2/TM4+TM2（Landst5-7） RVI=TM5-TM2/TM5+TM2（Landst8-9）GNDVI=NIR-G/NIR+GNIR是近红外波段，R是绿光波段的反射率。2.Python代码实现运行环境 arcpy python2.7 ILDE GUI运行。使用数据 Landsat8 OLI 未修正代码"""@soderayer绿度归一化植被指数GNDVI=TM4-TM2/TM4+TM2 G...

1.比值植被指数RVI公式：RVI=TM4/TM3（Landst5-7） RVI=TM5/TM4（Landst8-9）RVI=NIR/RNIR是近红外波段，R是红光波段的反射率。对于浓密植物反射的红光辐射很小,RVI将无限增长。该指数可以反映植被在可见光、近红外波段反射与土壤背景之间差异的指标，在一定条件下能用来定量说明植被的生长状况。1．绿色健康植被覆盖地区RVI远大于1，而无植被覆盖的地面（裸土、人工建筑、水体、植被枯死或严重虫害）的RVI在1附近。植被的RVI通常大于2；..

DEM(数字高程模型)数据可以计算不同的地形指数，地形指数也有多种，可以表达不同的物理意义。本文主要计算指定窗口大小下，地形高差指数，比如3*3窗口下，最大值和最小值之差，计算方法也比较简单。其他地形指数还有很多，比如地形耐用指数（Terrain Ruggedness Index，TPI）、地形位置指数（Topographic Position Index，TRI）和粗糙度（roughness).

1.叶面积指数概念 叶面积指数（leaf area index, LAI）不同资料显示定义有所不同，有的说是指单位地面上的绿叶面积，是植被冠层最显著的特征之一，具体是指一定土地面积上植物叶面面积总和与土地面积之比。它是植被光合作用模型和蒸散模型中的重要参数之一。有的说是是一块地上阳光直射时作物叶片垂直投影的总面积与占地面积的比值。即：叶面积指数=投影总面积/占地面积。2.叶面积指数计算方法同上一篇文章，使用遥感反演方法反演LAI,具体公式如下：公式：LAI=3.618*EVI-0.1..

1.叶面积指数概念 叶面积指数（leaf area index, LAI）不同资料显示定义有所不同，有的说是指单位地面上的绿叶面积，是植被冠层最显著的特征之一，具体是指一定土地面积上植物叶面面积总和与土地面积之比。它是植被光合作用模型和蒸散模型中的重要参数之一。有的说是是一块地上阳光直射时作物叶片垂直投影的总面积与占地面积的比值。即：叶面积指数=投影总面积/占地面积。2.叶面积指数计算方法 叶面积指数测定方法有多种，比如百度百科中将其分为直接方法和间接方法。直接测定方法是一种传统的、...

"""@soderayer date:20220428to read NDVI filestestfile Landsat8 collection2 level2 NDVI datacalculate the confidence interval to csvdata info table：Attribute ValueLong Name Normalized Difference Vegetation IndexShort Name LC8NDVI, LE7NDVI, L...

1.MOD09A1数据介绍1.1 数据介绍网址：LP DAAC - Search中分辨率成像光谱仪 (MODIS) Terra MOD09A1 6.1 版产品提供了对 Terra MODIS 波段 1 至 7 的表面光谱反射率的估计，并针对大气条件（如气体、气溶胶和瑞利散射）进行了校正。 除了七个 500 米 (m) 反射带之外，还有两个质量层和四个观测带。 对于每个像素，从 8 天合成周期内的所有采集中选择一个值。 像素选择的标准包括云天顶和太阳天顶。 当多次采集满足标准时，使用具有最小通道 3（

1.MOD11A2数据介绍MOD11A2 第 6 版产品在 1,200 x 1,200 公里的网格中以 1 公里 (km) 的空间分辨率提供平均 8 天的每像素地表温度和发射率 (LST&E)。 MOD11A2 中的每个像素值是在该 8 天内收集的所有相应 MOD11A1 LST 像素的简单平均值。 选择 8 天的合成周期是因为该周期是 Terra 和 Aqua 平台的精确地面轨道重复周期的两倍。 与白天和夜间地表温度带一起提供的是相关的质量控制评估、观测时间、视角天顶角和晴空覆盖率以及来自土地

"""@soderayer date:20220428to read NDVI filestestfile Landsat8 collection2 level2 NDVI datacalculate the confidence interval to csvdata info table：Attribute ValueLong Name Normalized Difference Vegetation IndexShort Name LC8NDVI, LE7ND...

我们在处理大批量压缩包数据的时候，有时候只需要解压压缩包里面指定内容的数据，为了少占用磁盘空间，我们可以使用python中的tarfile库来批量解压大批量压缩包数据中指定数据到指定文件夹里面。1.使用的python库import tarfile2.使用数据下载的Landsat8压缩包数据,文件名称如：LC081230382020041302T1-SC20220427023103.tar.gztf = tarfile.open(filename, "r:gz")filenames

1.数据下载在EarthExplorer查询网站提交订单后，发现一直无法下载，因此，将订单信息做成产品名称列表，将查询结果导出为csv文件格式，然后将产品名称一列做成txt文本列表，（如下图，导出元数据）在ESPA - LSRD提交这个文本订单并下载，具体过程网上有不少。点击New Order菜单进入提交产品列表页面，这里可以将数据转换为自己需要的投影和分辨率。更多处理选项里面还可以下载反射率产品和NDVI等指数产品，提交的时候报错了，就只下载了原始数据产品。提交后就开始处理..