基于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\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, 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\n", *longtitude); //printf("中央纬度为:%f\n", *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\n",*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\n", *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\n", *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\n", *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\n", *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\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()[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开发~~~~
