基于Google Earth Engine的全国地表温度反演

国内研究landsat8温度反演的人员很多，但是现有算法一般都是一景为例子，进行开展。

这有一个局限性，当研究的尺度很大时，就需要比较大的运算量了，例如全省温度，全国温度，全球温度，当然大家可能会说，可以用modis数据进行，我们当然不用modis，我用Landsat做中高尺度温度反演，我目前基于单通道算法，做出了全国中高分辨率的温度产品，目前来说，效果还是不错的，跟大气传输方法进行了对比，精度控制在1摄氏度以内。

我选择2017年进行实验，GEE总共运算了4759景Landsat影像，耗时大概1.3s，对比与传统的本地数据处理，一个天上一个地下，在这里，再次强烈建议我国也应当搭建遥感大数据云计算平台，把我国的高分1号，高分2号等影像数据，都集成到自己平台，然后提供API给相关人员调用，这样也可以免费，也可以收费，商业服务也上了一个档次。

不多说，我们看一下成果：

总体来看视觉效果真的很震撼，我计算的全年均值结果，可以看到新疆除了天山区域，大部分地区都是高温区域，而西藏区域常年低温，宁夏地区也是常年温度较高。

看到上面的结果比较理想，应该可以进行工程化应用了，这时候当然得上我们坑王-c++编程来实现工程化了。。。

此处省略一万字，我们重点参考这篇博文--http://blog.sina.com.cn/s/blog_764b1e9d0102wbl7.html，在此感谢邓书斌老师的ENVI教程，学习良多，不多，贴出我们的c++，opensourCode

第一段代码定义我们的头文件，方便后面调用：

//#include "stdafx.h"

#include "math.h"

#include <iostream>

#include <fstream>

#include <cassert>

#include <string>

#include <sstream>

//#include "gdal.h"

//#include "gdal_priv.h"

#include <cstdlib>

using namespace std;

// Landsat8温度反演头文件

class Landsat8_LST {

public:

    struct imgData

    {

        double adfGeoTransform[];  // 变换参数

        float *pData; // 影像数据,统一转换为float类型方便像元处理

        int imglength;// 图像大小

        int imgWidth;

        int imgHeight;

        int imgBandNum;

    };

    // 辐射定标

    void radImg(float *band, float add, float bias, int imglength) ;

    // 黑体反射率

    float* SurEmiss(float *nir, float *red, int imglength) ;

    // 温度反演

    float* LST(float t, float Lu, float Ld, float *aima, float *band10_rad, int imglength) ;

    // 读取影像

    void readImage(char *imgpath, imgData *IMG) ;

    // 写出影像

    void writeImage(char *imgPath, float *Img, int nImgSizeX, int  nImgSizeY, int nBandCount, double *adfGeoTransform) ;

    //去掉字符串中的空值

    string removeNUllString(string sNewTag) ;

    //string转换为其他类型

    template <class Type>

    Type stringToNum(const string &str) ;

    //读取Landsat8源文件

    void getMetaData(string metaDataPath, float *R_add, float *NIR_add, float *TIRS_add,

        float *R_mult, float *NIR_mult, float *TIRS_mult, string getTime,

        float *longtitude, float *latitude) ;

    //Landsat8温度反演

    void LandSat8Image_temperatureRetreval(char *RedBandPath, char *NirBandPath, char *TIRSBandPath,

        float t, float Lu, float Ld, string metaDataPath, char *imgPath) ;

};

接着是类的具体实现：

#include "stdafx.h"

#include "math.h"

#include <iostream>

#include <fstream>

#include <cassert>

#include <string>

#include <sstream>

#include "gdal.h"

#include "gdal_priv.h"

#include <cstdlib>

#include "LandSurfaceTemp.h"

using namespace std;

/*Landsat8 波段辐射定标

* add: 增益参数 bias:偏移参数

*/

void Landsat8_LST::radImg(float *band, float add, float bias, int imglength) {

    for (int i = ; i < imglength; i++) {

        band[i] = band[i] * bias - add;

    }

}

/*比辐射率计算

* nir: 近红外波段 red:红波段

*/

float*  Landsat8_LST::SurEmiss(float *nir, float *red, int imglength) {

    float *FV = new float[imglength];

    float *NDVI = new float[imglength];

    float *aima = new float[imglength];

    for (int i = ; i < imglength; i++)

    {

        NDVI[i] = (nir[i] - red[i]) / (nir[i] + red[i]);

        if (NDVI[i] > 0.7)

        {

            FV[i] = ;

        }

        if (NDVI[i]<)

        {

            FV[i] = ;

        }

        if ((NDVI[i] < 0.7) && (NDVI[i] > ))

        {

            FV[i] = (NDVI[i] - 0.05) / (0.7 - 0.05);

        }

        aima[i] = FV[i] * 0.004 + 0.986;

    }

    delete[] FV;

    delete[] NDVI;

    return aima;

}

/*地表温度反演

* t: 大气参数

* Lu: 大气参数

* Ld: 大气参数

* aima: 地表比辐射率

* band10_rad: 辐射定标后第10波段

*/

float* Landsat8_LST::LST(float t, float Lu, float Ld, float *aima, float *band10_rad, int imglength) {

    float *LST = new float[imglength];

    float *BlaRad = new float[imglength];

    for (int i = ; i < imglength; i++)

    {

        BlaRad[i] = (band10_rad[i] - Lu - t*( - aima[i]) * Ld) / (t*aima[i]);

        LST[i] = 1321.08 / log(774.89 / BlaRad[i] + ) - ;

    }

    delete[] BlaRad;  //释放内存

    delete[] aima;

    return LST;

}

/*栅格影像读取，返回数据指针

* imgPath:图像位置

* 返回float类型的数据指针

*/

void Landsat8_LST::readImage(char *imgpath, imgData *IMG) {

    GDALDataset *img = (GDALDataset*)GDALOpen(imgpath, GA_ReadOnly);

    if (img != NULL) {

        int imgWidth = img->GetRasterXSize();   //图像宽度，特别注意：对应matlab中的行

        int imgHeight = img->GetRasterYSize();  //图像高度，特别注意：对应matlab中的列

        int bandNum = img->GetRasterCount();    //波段数

        IMG->imgBandNum = bandNum;

        IMG->imgHeight = imgHeight;

        IMG->imgWidth = imgWidth;

        GDALRasterBand *poBand;

        poBand = img->GetRasterBand();  //灰度一个波段

                                         //printf("影像数据类型为:%s\n", GDALGetDataTypeName(poBand->GetRasterDataType()));

        img->GetGeoTransform(IMG->adfGeoTransform); // 变换参数

        int size = imgWidth*imgHeight*bandNum;

        IMG->pData = new float[size]; //分配缓冲区空间

        IMG->imglength = size;

        //读取

        poBand->RasterIO(GF_Read, , , imgWidth, imgHeight, IMG->pData,

            imgWidth, imgHeight, GDT_Float32, , );

        GDALClose(img); // 释放内存

    }

}

/*写出栅格影像

* imgPath:输出影像位置

* adfGeoTransform:变换参数

* IMG:导出的影像数组

*/

void Landsat8_LST::writeImage(char *imgPath, float *Img, int nImgSizeX, int  nImgSizeY, int nBandCount, double *adfGeoTransform) {

    GDALDataset *poDataset2; //待创建的GDAL数据集

    GDALDriver *poDriver;    //驱动，用于创建新的文件

                             //创建新文件

    poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

    //获取格式类型

    char **papszMetadata = poDriver->GetMetadata();

    //特别注意，数据类型要与后面的写出类型要保持一致

    poDataset2 = poDriver->Create(imgPath, nImgSizeX, nImgSizeY, nBandCount, GDT_Float32, papszMetadata);

    //坐标赋值

    poDataset2->SetGeoTransform(adfGeoTransform);

    //将图像数据写入新图像中

    poDataset2->RasterIO(GF_Write, , , nImgSizeX, nImgSizeY,

        Img, nImgSizeX, nImgSizeY, GDT_Float32, nBandCount, , , , );

    GDALClose(poDataset2);

    delete poDriver;

}

/*读取Landsat8元数据中的辐射定标参数，拍摄时间，中央经纬度信息

* metaDataPath:元数据路径

* R_add: 红波段增益 ;

* NIR_add 近红外波段增益;

* TIRS_add 第10波段增益;

* R_mult 红波段偏移;

* NIR_mult 近红外波段偏移;

* TIRS_mult 第10波段偏移;

* getTime:影像拍摄时间;

* longtitude:中央经度;

* latitude：中央纬度;

*/

// 去掉一行字符串中的空字符

string Landsat8_LST::removeNUllString(string sNewTag) {

    int begin = ;

    begin = sNewTag.find(" ", begin);

    while (begin != -)

    {

        sNewTag.replace(begin, , "");

        begin = sNewTag.find(" ", begin);

    }

    return sNewTag;

}

//模板函数：将string类型变量转换为常用的数值类型（此方法具有普遍适用性）

template <class Type>

Type Landsat8_LST::stringToNum(const string &str)

{

    istringstream iss(str);

    Type num;

    iss >> num;

    return num;

}

//读取Landsat8影像元数据

void Landsat8_LST::getMetaData(string metaDataPath, float *R_add, float *NIR_add, float *TIRS_add,

    float *R_mult, float *NIR_mult, float *TIRS_mult, string getTime,

    float *longtitude, float *latitude) {

    std::ifstream infile;

    infile.open(metaDataPath.data());   //将文件流对象与文件连接起来

    assert(infile.is_open());   //若失败,则输出错误消息,并终止程序运行 

    string s; //存储每一行的文本数据

    float ur_LAT = , ur_LON = , ll_LAR = , ll_LON = ; //右上和左下经纬度

    while (getline(infile, s))

    {

        string S = removeNUllString(s);

        // 得到拍摄时间

        string sen = "DATE_ACQUIRED=";

        int time = S.find(sen);

        if (time != -) {

            getTime = S.substr(sen.length(), S.length());

        }

        // 根据四个角度的经纬度，计算中央经纬度

        string UR_LAT = "CORNER_UR_LAT_PRODUCT=";

        time = S.find(UR_LAT);

        if (time != -) {

            string temp = S.substr(UR_LAT.length(), S.length());

            ur_LAT = stringToNum<float>(temp);

        }

        string UR_LON = "CORNER_UR_LON_PRODUCT=";

        time = S.find(UR_LON);

        if (time != -) {

            string temp = S.substr(UR_LON.length(), S.length());

            ur_LON = stringToNum<float>(temp);

        }

        string LL_LAR = "CORNER_LL_LAT_PRODUCT=";

        time = S.find(LL_LAR);

        if (time != -) {

            string temp = S.substr(LL_LAR.length(), S.length());

            ll_LAR = stringToNum<float>(temp);

        }

        string LL_LON = "CORNER_LL_LON_PRODUCT=";

        time = S.find(LL_LON);

        if (time != -) {

            string temp = S.substr(LL_LON.length(), S.length());

            ll_LON = stringToNum<float>(temp);

        }

        if ((ll_LAR != ) && (ur_LAT != ) && (ll_LON != ) && (ur_LON != )) {

            *latitude = (ll_LAR + ur_LAT) / ;

            *longtitude = (ll_LON + ur_LON) / ;

            //printf("中央经度为:%f\n", *longtitude);

            //printf("中央纬度为:%f\n", *latitude);

        }

        // 得到 R波段辐射定标增益参数

        string senR = "RADIANCE_ADD_BAND_4=";

        time = S.find(senR);

        if (time != -) {

            string temp = S.substr(senR.length(), S.length());

            *R_add = stringToNum<float>(temp);

            //printf("R波段辐射定标增益参数:%f\n",*R_add);

        }

        // 得到 NIR波段辐射定标增益参数

        string senNIR = "RADIANCE_ADD_BAND_5=";

        time = S.find(senNIR);

        if (time != -) {

            string temp = S.substr(senNIR.length(), S.length());

            *NIR_add = stringToNum<float>(temp);

            //printf("NIR波段辐射定标增益参数:%f\n", *NIR_add);

        }

        // 得到 TIRS波段辐射定标增益参数

        string senTIRS = "RADIANCE_ADD_BAND_10=";

        time = S.find(senTIRS);

        if (time != -) {

            string temp = S.substr(senTIRS.length(), S.length());

            *TIRS_add = stringToNum<float>(temp);

            //printf("TIRS波段辐射定标增益参数:%f\n", *TIRS_add);

        }

        // 得到 R波段辐射定标偏移参数

        string senRmult = "RADIANCE_MULT_BAND_4=";

        time = S.find(senRmult);

        if (time != -) {

            string temp = S.substr(senRmult.length(), S.length());

            *R_mult = stringToNum<float>(temp);

            //printf("R波段辐射定标偏移参数:%f\n", *R_mult);

        }

        // 得到 NIR波段辐射定标偏移参数

        string senNIRmult = "RADIANCE_MULT_BAND_5=";

        time = S.find(senNIRmult);

        if (time != -) {

            string temp = S.substr(senNIRmult.length(), S.length());

            *NIR_mult = stringToNum<float>(temp);

            //printf("NIR波段辐射定标偏移参数:%f\n", *NIR_mult);

        }

        // 得到 TIRS波段辐射定标偏移参数

        string senTIRSmult = "RADIANCE_MULT_BAND_10=";

        time = S.find(senTIRSmult);

        if (time != -) {

            string temp = S.substr(senTIRSmult.length(), S.length());

            *TIRS_mult = stringToNum<float>(temp);

            //printf("TIRS波段辐射定标偏移参数:%f\n", *TIRS_mult);

        }

    }

    infile.close();             //关闭文件输入流

}

//基于大气校正法的Landat8地表温度反演

void Landsat8_LST::LandSat8Image_temperatureRetreval(char *RedBandPath, char *NirBandPath, char *TIRSBandPath,

    float t, float Lu, float Ld, string metaDataPath, char *imgPath) {

    // 读取红波段

    struct imgData *redband = new imgData;

    readImage(RedBandPath, redband);

    // 读取近红波段

    struct imgData *nirband = new imgData;

    readImage(NirBandPath, nirband);

    // 读取第10波段

    struct imgData *TIRSband = new imgData;

    readImage(TIRSBandPath, TIRSband);

    // 读取影像元数据参数

    string getTime;

    float *R_add = new float;  //记得初始化局部变量指针

    float *NIR_add = new float;

    float *TIRS_add = new float;

    float *R_mult = new float;

    float *NIR_mult = new float;

    float *TIRS_mult = new float;

    float *longtitude = new float;

    float *latitude = new float;

    getMetaData(metaDataPath, R_add, NIR_add, TIRS_add, R_mult, NIR_mult,

        TIRS_mult, getTime, longtitude, latitude);

    // 辐射定标

    radImg((float *)(redband->pData), *R_add, *R_mult, redband->imglength); //R波段辐射定标

    radImg((float *)(nirband->pData), *NIR_add, *NIR_mult, redband->imglength); //NIR波段辐射定标

    radImg((float *)(TIRSband->pData), *TIRS_add, *TIRS_mult, redband->imglength); //热红外波段辐射定标

                                                                                   // 比辐射率计算

    float *aima = SurEmiss(nirband->pData, redband->pData, redband->imglength);

    // 温度反演

    float *LandSurfaceTemp = LST(t, Lu, Ld, aima, TIRSband->pData, redband->imglength);

    //导出结果

    writeImage(imgPath, LandSurfaceTemp, redband->imgWidth,

        redband->imgHeight, redband->imgBandNum, redband->adfGeoTransform);

    // 清空内存

    delete[] redband->pData; //注意,使用了new数组，需要进行内存释放，对于结构体指针，先释放结构中的数组指针，然后释放该结构体

    delete[] nirband->pData;

    delete[] TIRSband->pData;

    delete redband, nirband, TIRSband;

    delete R_add, NIR_add, TIRS_add, R_mult, NIR_mult, TIRS_mult, latitude;

}

然后是我们的qt主函数实现：

#include "LandSurfaceTemp.h"

#include "QtGuiApplication.h"

#include "QMessageBox"

#include "QFileDialog"

#include <string>

#include "math.h"

#include <iostream>

#include <fstream>

#include <cassert>

#include <string>

#include <sstream>

#include "gdal.h"

#include "gdal_priv.h"

#include <cstdlib>

#pragma execution_character_set("utf-8") // 添加此句，支持中文

using namespace std;

//解决Qstring转string中文乱码函数

QString str2qstr(string str)

{

    return QString::fromLocal8Bit(str.data());

}

string qstr2str(QString qstr)

{

    QByteArray cdata = qstr.toLocal8Bit();

    return string(cdata);

}

//模板函数：将string类型变量转换为常用的数值类型（此方法具有普遍适用性）

template <class Type>

Type stringToNum(const string &str)

{

    istringstream iss(str);

    Type num;

    iss >> num;

    return num;

}

QtGuiApplication::QtGuiApplication(QWidget *parent)

    : QMainWindow(parent)

{

    ui.setupUi(this);

}

void QtGuiApplication::Open_Rband_buttonClick()  //打开红波段影像

{

    QFileDialog *fileDialog = new QFileDialog(this);

    fileDialog->setWindowTitle(tr("打开红波段影像"));

    fileDialog->setDirectory(".");

    if (fileDialog->exec() == QDialog::Accepted)

    {

        QString path = fileDialog->selectedFiles()[];

        ui.openImgEdit->setText(path); // 所有控件全都通过ui来调用

    }

    else

    {

        QMessageBox::information(NULL, tr("Path"), tr("您没有选择任何文件!"));

    }

}

void QtGuiApplication::Open_NIRband_buttonClick()  //打开近红波段影像

{

    QFileDialog *fileDialog = new QFileDialog(this);

    fileDialog->setWindowTitle(tr("打开近红波段影像"));

    fileDialog->setDirectory(".");

    if (fileDialog->exec() == QDialog::Accepted)

    {

        QString path = fileDialog->selectedFiles()[];

        ui.openImgEdit_2->setText(path); // 所有控件全都通过ui来调用

    }

    else

    {

        QMessageBox::information(NULL, tr("Path"), tr("您没有选择任何文件!"));

    }

}

void QtGuiApplication::Open_band10_buttonClick()  //打开第10波段影像

{

    QFileDialog *fileDialog = new QFileDialog(this);

    fileDialog->setWindowTitle(tr("打第10波段影像"));

    fileDialog->setDirectory(".");

    if (fileDialog->exec() == QDialog::Accepted)

    {

        QString path = fileDialog->selectedFiles()[];

        ui.openImgEdit_3->setText(path); // 所有控件全都通过ui来调用

    }

    else

    {

        QMessageBox::information(NULL, tr("Path"), tr("您没有选择任何文件!"));

    }

}

void QtGuiApplication::OpenMetaFileClick()  //打开源文件

{

    QFileDialog *fileDialog = new QFileDialog(this);

    fileDialog->setWindowTitle(tr("打开源文件"));

    fileDialog->setDirectory(".");

    if (fileDialog->exec() == QDialog::Accepted)

    {

        QString path = fileDialog->selectedFiles()[];

        ui.openImgEdit_4->setText(path); // 所有控件全都通过ui来调用

    }

    else

    {

        QMessageBox::information(NULL, tr("Path"), tr("您没有选择任何文件!"));

    }

}

void QtGuiApplication::savebuttonclick()//保存温度反演影像

{

    QString fileName;

    fileName = QFileDialog::getSaveFileName(this,

        tr("Save Image!"), "", tr("Img File (*.tif)"));

    if (!fileName.isNull())

    {

        ui.saveImgEdit->setText(fileName);

    }

    else

    {

        QMessageBox::information(NULL, tr("Path"), tr("您没有选择任何文件!"));

    }

}

//执行按钮

void QtGuiApplication::LSTButtonClick()

{

    try

    {

        //必须先注册一个！

        GDALAllRegister();

        CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

        //打开红波段影像

        string Rband_imgPath = qstr2str(ui.openImgEdit->toPlainText());

        char *Rband_ImgPath = (char*)Rband_imgPath.data();

        //打开近红波段影像

        string NIRband_imgPath = qstr2str(ui.openImgEdit_2->toPlainText());

        char *NIRband_ImgPath = (char*)NIRband_imgPath.data();

        //打开第10波段影像

        string band10_imgPath = qstr2str(ui.openImgEdit_3->toPlainText());

        char *band10_ImgPath = (char*)band10_imgPath.data();

        //打开元数据

        string meta_imgPath = qstr2str(ui.openImgEdit_4->toPlainText());

        char *meta_ImgPath = (char*)meta_imgPath.data();

        //大气参数

        string t = ui.T_Edit->text().toStdString();

        string Ld = ui.Ld_Edit_2->text().toStdString();

        string Lu = ui.Lu_Edit_3->text().toStdString();

        float T = stringToNum<float>(t);

        float LD = stringToNum<float>(Ld);

        float LU = stringToNum<float>(Lu);

        //保存文件

        string sPath = qstr2str(ui.saveImgEdit->toPlainText());

        char *savePath = (char*)sPath.data();

        // 温度反演

        Landsat8_LST lst; //新建类

        lst.LandSat8Image_temperatureRetreval(Rband_ImgPath, NIRband_ImgPath, band10_ImgPath,

            T, LU, LD, meta_ImgPath, savePath);

        QMessageBox::information(NULL, tr("提示！"), tr(" 温度反演成功！"));

    }

    catch (const char* msg)

    {

        QMessageBox::information(NULL, tr("错误！"), tr("出现异常"));

    }

}

然后是画出我们的界面：

上篇文章写的全国水体研究，大家都比较喜欢，这篇全国温度反演，跟上篇属于姊妹篇，一直没有时间写，这次终于有时间写完了，如有代码和成果需求，请加qq1044625113，请注明，全国温度反演！

打个小广告：本人兼职GEE开发~~~~

基于Google Earth Engine的全国地表温度反演

基于Google Earth Engine的全国地表温度反演的相关教程结束。

相关推荐

[模板] 容斥原理: 二项式反演 / Stirling 反演 / min-max 容斥 / 子集反演 / 莫比乌斯反演

bzoj5093图的价值：多项式，斯特林数（二项式反演）

2022-12-09：上升的温度。以下的数据输出2和4，2015-01-02 的温度比前一天高（10 -＞ 25），2015-01-04 的温度比前一天高（20 -＞ 30），sql语句如何写？ DR

莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛学习笔记

单窗算法的地表温度反演：谷歌地球引擎GEE代码

Note -「单位根反演」学习笔记

2022年写的香橙派 OrangePi Zero 用python获取dht11温度和湿度

单窗算法的地表温度反演：谷歌地球引擎GEE代码

基于Google Earth Engine的全国地表温度反演

基于Google Earth Engine的全国地表温度反演的相关教程结束。

相关推荐

[模板] 容斥原理: 二项式反演 / Stirling 反演 / min-max 容斥 / 子集反演 / 莫比乌斯反演

bzoj5093图的价值：多项式，斯特林数（二项式反演）

2022-12-09：上升的温度。以下的数据输出2和4，2015-01-02 的温度比前一天高（10 -＞ 25），2015-01-04 的温度比前一天高（20 -＞ 30），sql语句如何写？ DR

莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛 学习笔记

单窗算法的地表温度反演：谷歌地球引擎GEE代码

Note -「单位根反演」学习笔记

2022年写的香橙派 OrangePi Zero 用python获取dht11温度和湿度

单窗算法的地表温度反演：谷歌地球引擎GEE代码

莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛学习笔记