Google Earth Engine城市水体提取

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("读取影像成功！
");

    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("提取水体成功！
");

    // 清空内存
    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）。