Google Earth Engine城市水体提取

2018/11/15 11:29
阅读数 472

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)。

展开阅读全文
打赏
0
0 收藏
分享
加载中
更多评论
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部