Java調用GDAL根據指定坐标點讀取栅格資料像元值

目錄

一、前言

二、實作思路 

三、實作的代碼

1、程式代碼

2、運作結果

3、驗證

四、補充說明

一、前言

在3S技術（全球衛星定位系統、地理資訊系統、遙感）中，栅格資料是一類常見的資料，栅格資料有其特定的資料結構。簡單點說，栅格資料的表現的資訊主要展現在每個像元對應的數值上，可以抽象為數學上的“矩陣”，這一點與數字圖像類似，與其不同的是栅格資料往往帶有空間資訊（地理坐标、投影坐标），是以根據坐标點讀取栅格資料的值是一項基礎的工作。

GDAL是開源的矢量、栅格資料的處理程式，支援較多的程式設計語言，但網絡上有關Java的文章相對C和Python來說很少。是以在這裡記錄一下，已備不時之需。

二、實作思路 

根據GDAL的API文檔，解析的思路應該是，首先根據檔案路徑擷取到對應的DataSet類，然後從DataSet類擷取Band類，然後調用Band類的ReadRaster擷取一個一維數組，這個一維數組就是該栅格資料的某一區域的值。當ReadRaster的參數範圍足夠小時，就可以擷取到一個像元的值，即一維數組隻有一個值。

由于ReadRaster的參數需要像素坐标來定位，是以需要先使用DataSet類的GetGeoTransform方法擷取像素坐标轉地圖坐标的參數，然後根據像素坐标轉地圖坐标的公式和指定的地圖坐标點，解二進制一次方程組，得到像素坐标。進而調用ReadRaster方法。

三、實作的代碼

1、程式代碼

package com.myself.raster;

import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;

import java.util.Arrays;

public class RasterDemo {

    public static void main(String[] args) {
        String filepath = "F:\\raster\\Demo.tif";
        //注冊
        gdal.AllRegister();
        //打開檔案擷取資料集
        Dataset dataset = gdal.Open(filepath,
                gdalconstConstants.GA_ReadOnly);
        if (dataset == null) {
            System.out.println("打開"+filepath+"失敗"+gdal.GetLastErrorMsg());
            System.exit(1);
        }
        //擷取驅動
        Driver driver = dataset.GetDriver();
        //擷取驅動資訊
        System.out.println("driver long name: " + driver.getLongName());
        //擷取栅格數量
        int bandCount = dataset.getRasterCount();
        System.out.println("RasterCount: " + bandCount);
        //構造仿射變換參數數組，并擷取資料
        double[] gt = new double[6];
        dataset.GetGeoTransform(gt);
        System.out.println("仿射變換參數"+Arrays.toString(gt));

        //指定經緯度
        double Latitude  = 86.053;
        double longitude = 16.529;

        //經緯度轉換為栅格像素坐标
        int[] ColRow=Coordinates2ColRow(gt,Latitude,longitude);

        //判斷是否坐标超出範圍
        if(ColRow[0]<0||ColRow[1]<0||ColRow[0]>dataset.getRasterXSize()||ColRow[1]>dataset.getRasterYSize()){
            System.out.println(Arrays.toString(ColRow)+"坐标值超出栅格範圍！");
            return;
        }

        //周遊波段，擷取該點對應的每個波段的值并列印到螢幕
        for (int i = 0;i<bandCount;i++){
            Band band = dataset.GetRasterBand(i+1);
            double[] values = new double[1];
            band.ReadRaster(ColRow[0], ColRow[1], 1, 1, values);
            double value = values[0];
            System.out.println("橫坐标："+Latitude +","+"縱坐标:"+longitude);
            System.out.format("Band"+(i+1)+": %s", value);
        }
        //釋放資源
        dataset.delete();

    }

    /**
     * 将地圖坐标轉換為栅格像素坐标
     * @param gt 仿射變換參數
     * @param X 橫坐标
     * @param Y 縱坐标
     * @return
     */
    public static int[] Coordinates2ColRow(double[] gt, double X, double Y){
        int[] ints = new int[2];
//        X = gt[0] + Xpixel*gt[1] + Yline*gt[2];
//        Y = gt[3] + Xpixel*gt[4] + Yline*gt[5];
//        消元法解二進制一次方程組
//        X-gt[0]=Xpixel*gt[1] + Yline*gt[2]
//        Xpixel = (X-gt[0] - Yline*gt[2])/gt[1]
//        Y - gt[3] = ((X-gt[0] - Yline*gt[2])/gt[1])gt[4] + Yline*gt[5]
//        (Y - gt[3])*gt[1] = ((X-gt[0] - Yline*gt[2]))*gt[4] + Yline*gt[5]*gt[1]
//        (Y - gt[3])*gt[1] =(X-gt[0])*gt[4] - Yline*gt[2]*gt[4] + Yline*gt[5]*gt[1]
//        (Y - gt[3])*gt[1] - (X-gt[0])*gt[4] = Yline(gt[5]*gt[1]-gt[2]*gt[4])

        //向下取整,如果向上取整會導緻計算結果偏大，進而在後面讀取到鄰近像元的資料
        double  Yline = Math.floor(((Y - gt[3])*gt[1] - (X-gt[0])*gt[4]) / (gt[5]*gt[1]-gt[2]*gt[4]));
        double  Xpixel = Math.floor((X-gt[0] - Yline*gt[2])/gt[1]);
        ints[0] = new Double(Xpixel).intValue();
        ints[1] = new Double(Yline).intValue();
        return ints;
    }

}
           

2、 運作結果

3、驗證

這裡使用QGIS來打開資料檔案，并使用識别功能來驗證有關點的值。

 可以發現，程式解析結果與QGIS解析結果一緻。

四、補充說明

dataset.GetGeoTransform擷取的gt[]含義如下：

//gt[0]  左上角x坐标

//gt[1]  東西方向分辨率

//gt[2]  旋轉角度1, 0表示圖像 "北方朝上"

//gt[3]  左上角y坐标

//gt[4]  旋轉角度2, 0表示圖像 "北方朝上"

//gt[5]  東西方向分辨率

地圖坐标X =  左上角x坐标 + 列數*東西方向分辨率 + 行數*旋轉角度1；

地圖坐标Y =  左上角y坐标 + 列數*旋轉角度2 + 行數*南北方向分辨率；

band.ReadRaster的入參含義如下：public int ReadRaster(int xoff, int yoff, int xsize, int ysize, double[] array)

public int ReadRaster(int xoff, int yoff, int xsize, int ysize, double[] array)

xoff 想要讀取的部分原點位置橫圖像坐标

yoff 想要讀取的部分原點位置縱圖像坐标

xsize 指定要讀取部分圖像的矩形邊長。 

ysize 指定要讀取部分圖像的矩形邊長。 

array 用來存儲讀取到的資料的數組。

本文的資料是作者自己做的，沒有任何業務意義。如有需要可以下載下傳：Demo.tif