Google Earth Engine城市水体提取
大家都知道城市水体提取相比较于山区，丘陵的地区，肯定是比较难的，为什么呢，因为城市水体有很多高层建筑导致的阴影，这个就非常复杂了，而且现在很多高分影像只有可见光和近红外波段，那么我们如何准确提取城市水体呢？
Remoe Sensing2018年刊发了一篇城市水体高分影像自动提取算法（Two-Step Urban Water Index (TSUWI): A New Technique for High-Resolution Mapping of Urban SurfaceWater [J]Remote Sensing，2018），初步看来，效果还行，在高分二号上面效果不错，我再想，如果对于开源的哨兵、Landsat如何？这些是中等分辨率影像，能做到吗？
话不多说，利用GEE，直接编码，实验结果如下（以2018年10月的北京某景Sentinel2影像为例）：
(a) 这是原始影像
(b) 这是城市水体指数
(c) 这是城市阴影指数
(d) 这是城市水体提取结果，蓝色为水体
其中城市水体指数和城市阴影指数计算公式如下所示：
我把最终成果发布成了APPengine(https://wang749195.users.earthengine.app/view/urbanwaterextraction)，大家可以直接在web上看，总的来说，实验结果还是不错的，去掉了阴影现象，这篇文章出自中科院遥感所，在此申明，值得一读，后续我会发布C++软件版本，Matlab版本，以及Python版本。我个人的开发思路是，首先用GEE实现，如果GEE不好实现，就用matlab或者python实现第一遍，效果可以，能工程应用，立马就用GDAL+C++打包成工程源代码，我感觉这样会节省时间，且不会造成时间浪费。
接着上面讲，我们用c++来实现一遍，使用GDAL读写影像，先把这两个函数写上来：
/*栅格影像读取，返回数据指针
* imgPath:图像位置
* 返回float类型的数据指针
*/
void readImage(char *imgpath, imgData *IMG, int bandindex) {
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->imgH = imgHeight;
IMG->imgW = imgWidth;
GDALRasterBand *poBand;
poBand = img->GetRasterBand(bandindex); //灰度一个波段
img->GetGeoTransform(IMG->adfGeoTransform); // 变换参数
int size = imgWidth*imgHeight;
IMG->pData = new float[size]; //分配缓冲区空间
//读取
poBand->RasterIO(GF_Read, 0, 0, imgWidth, imgHeight, IMG->pData,
imgWidth, imgHeight, GDT_Float32, 0, 0);
GDALClose(img); // 释放内存
}
}
/*写出栅格影像
* imgPath:输出影像位置
* adfGeoTransform:变换参数
* IMG:导出的影像数组
*/
void 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;
}然后就是我们的USI,UWI计算公式，贴上来：
// 计算UWI指数
void UWI_cal(float *rband, float *gband, float *nirband,float *UWI,int width,int length) {
int Length = width*length;
for (int i = 0; i < Length; i++) {
UWI[i] = (gband[i] - 1.1*rband[i] - 5.2*nirband[i] + 0.4) /
abs(gband[i] - 1.1*rband[i] - 5.2*(nirband[i]));
}
}
// 计算USI指数
void USI_cal(float *rband, float *gband, float *bband, float *nirband, float *USI, int width, int length) {
int Length = width*length;
for (int i = 0; i < Length; i++) {
USI[i] = 0.25*gband[i] / rband[i] - 0.57*nirband[i] /
gband[i] - 0.83*bband[i] / gband[i] + 1.0;
}
}然后就是我们的影像数据结构：
/*可见光与近红外波段数据结构
*/
struct imgData {
float *pData;
int imgH;
int imgW;
double adfGeoTransform[6];
};last but not least，就是我们的main函数：
int main()
{
//必须先注册一个！
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
char *ImgPath = "C:\\Users\\Administrator\\Desktop\\UrbanWater\\SentinelImg.tif";
// 读取蓝波段
imgData *B = new imgData;
readImage(ImgPath, B, 1);
// 读取绿波段
imgData *G = new imgData;
readImage(ImgPath, G, 2);
// 读取红波段
imgData *R = new imgData;
readImage(ImgPath, R, 3);
// 读取近红外波段
imgData *NIR = new imgData;
readImage(ImgPath, NIR, 4);
printf("读取影像成功！\n");
int width = B->imgW;
int height = B->imgH;
float *USI = new float[width*height];
float *UWI = new float[width*height];
UWI_cal(R->pData, G->pData, NIR->pData, UWI, width, height);
USI_cal(R->pData, G->pData, B->pData, NIR->pData, USI, width, height);
float T1 = -0.1;
float T2 = -0.2;
float *UrbanWater = new float[width*height];
UrbanWaterExtraction(T1, T2, UWI, USI, UrbanWater, width, height);
char *savePath = "C:\\Users\\Administrator\\Desktop\\UrbanWater\\urbanwater.tif";
writeImage(savePath, UrbanWater, width, height, 1, R->adfGeoTransform);
printf("提取水体成功！\n");
// 清空内存
delete []NIR->pData;
delete []R->pData;
delete []G->pData;
delete []B->pData;
delete []UrbanWater;
delete []USI;
delete []UWI;
delete NIR, R, G, B;
system("pause");
}还是上一张c++搞出来的城市水体图吧：
 |
可以看到，GEE与c++效果几乎一样，但是GEE的栅格渲染，还是非常值得国产软件学习！
（打个小广告，本文兼职软件开发，qq1044625113）。
