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