基于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[6];  // 变换参数
        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 = 0; 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 = 0; i < imglength; i++)
    {
        NDVI[i] = (nir[i] - red[i]) / (nir[i] + red[i]);
        if (NDVI[i] > 0.7)
        {
            FV[i] = 1;
        }
        if (NDVI[i]<0)
        {
            FV[i] = 0;
        }
        if ((NDVI[i] < 0.7) && (NDVI[i] > 0))
        {
            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 = 0; i < imglength; i++)
    {
        BlaRad[i] = (band10_rad[i] - Lu - t*(1 - aima[i]) * Ld) / (t*aima[i]);
        LST[i] = 1321.08 / log(774.89 / BlaRad[i] + 1) - 273;
    }
    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(1);  //灰度一个波段
                                         //printf("影像数据类型为:%s
", GDALGetDataTypeName(poBand->GetRasterDataType()));

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

        int size = imgWidth*imgHeight*bandNum;
        IMG->pData = new float[size]; //分配缓冲区空间
        IMG->imglength = size;

        //读取
        poBand->RasterIO(GF_Read, 0, 0, imgWidth, imgHeight, IMG->pData,
            imgWidth, imgHeight, GDT_Float32, 0, 0);

        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, 0, 0, nImgSizeX, nImgSizeY,
        Img, nImgSizeX, nImgSizeY, GDT_Float32, nBandCount, 0, 0, 0, 0);

    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 = 0;

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

    while (begin != -1)
    {

        sNewTag.replace(begin, 1, "");

        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 = 0, ur_LON = 0, ll_LAR = 0, ll_LON = 0; //右上和左下经纬度
    while (getline(infile, s))
    {
        string S = removeNUllString(s);

        // 得到拍摄时间
        string sen = "DATE_ACQUIRED=";
        int time = S.find(sen);
        if (time != -1) {
            getTime = S.substr(sen.length(), S.length());
        }

        // 根据四个角度的经纬度，计算中央经纬度
        string UR_LAT = "CORNER_UR_LAT_PRODUCT=";
        time = S.find(UR_LAT);
        if (time != -1) {
            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 != -1) {
            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 != -1) {
            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 != -1) {
            string temp = S.substr(LL_LON.length(), S.length());
            ll_LON = stringToNum<float>(temp);
        }

        if ((ll_LAR != 0) && (ur_LAT != 0) && (ll_LON != 0) && (ur_LON != 0)) {
            *latitude = (ll_LAR + ur_LAT) / 2;
            *longtitude = (ll_LON + ur_LON) / 2;
            //printf("中央经度为:%f
", *longtitude);
            //printf("中央纬度为:%f
", *latitude);
        }

        // 得到 R波段辐射定标增益参数
        string senR = "RADIANCE_ADD_BAND_4=";
        time = S.find(senR);
        if (time != -1) {
            string temp = S.substr(senR.length(), S.length());
            *R_add = stringToNum<float>(temp);
            //printf("R波段辐射定标增益参数:%f
",*R_add);
        }

        // 得到 NIR波段辐射定标增益参数
        string senNIR = "RADIANCE_ADD_BAND_5=";
        time = S.find(senNIR);
        if (time != -1) {
            string temp = S.substr(senNIR.length(), S.length());
            *NIR_add = stringToNum<float>(temp);
            //printf("NIR波段辐射定标增益参数:%f
", *NIR_add);
        }

        // 得到 TIRS波段辐射定标增益参数
        string senTIRS = "RADIANCE_ADD_BAND_10=";
        time = S.find(senTIRS);
        if (time != -1) {
            string temp = S.substr(senTIRS.length(), S.length());
            *TIRS_add = stringToNum<float>(temp);
            //printf("TIRS波段辐射定标增益参数:%f
", *TIRS_add);
        }

        // 得到 R波段辐射定标偏移参数
        string senRmult = "RADIANCE_MULT_BAND_4=";
        time = S.find(senRmult);
        if (time != -1) {
            string temp = S.substr(senRmult.length(), S.length());
            *R_mult = stringToNum<float>(temp);
            //printf("R波段辐射定标偏移参数:%f
", *R_mult);
        }

        // 得到 NIR波段辐射定标偏移参数
        string senNIRmult = "RADIANCE_MULT_BAND_5=";
        time = S.find(senNIRmult);
        if (time != -1) {
            string temp = S.substr(senNIRmult.length(), S.length());
            *NIR_mult = stringToNum<float>(temp);
            //printf("NIR波段辐射定标偏移参数:%f
", *NIR_mult);
        }

        // 得到 TIRS波段辐射定标偏移参数
        string senTIRSmult = "RADIANCE_MULT_BAND_10=";
        time = S.find(senTIRSmult);
        if (time != -1) {
            string temp = S.substr(senTIRSmult.length(), S.length());
            *TIRS_mult = stringToNum<float>(temp);
            //printf("TIRS波段辐射定标偏移参数:%f
", *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()[0];
        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()[0];
        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()[0];
        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()[0];
        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开发~~~~