Java calls GDAL to read raster data pixel values according to the specified coordinate points

Posted by fallen_angel21 on Tue, 05 Oct 2021 22:08:52 +0200

catalogue

1, Foreword

2, Realization idea  

3, Implemented code

1. Program code

2. Operation results

3. Verify

4, Supplementary notes

1, Foreword

In 3S Technology (global positioning system, geographic information system, remote sensing), grid data is a kind of common data, and grid data has its specific data structure. In short, the information represented by grid data is mainly reflected in the corresponding value of each pixel and can be abstracted into a mathematical "matrix", which is similar to digital images. The difference is that grid data often contains spatial information (geographic coordinates and projection coordinates), so reading the value of grid data according to coordinate points is a basic work.

GDAL is an open source vector and raster data processing program, which supports many programming languages, but there are few posts about Java on the network compared with C and Python. So record it here. It's ready for a rainy day.

2, Realization idea  

according to API documentation for GDAL The analytic way of thinking is to first get the corresponding DataSet class according to the file path, then get the Band class from the DataSet class, and then call the Band class ReadRaster to get a one-dimensional array, which is the value of a certain area of the raster data. When the parameter range of readraster is small enough, you can obtain the value of a pixel, that is, the one-dimensional array has only one value.

Because the parameters of ReadRaster need pixel coordinates to locate, you need to use the GetGeoTransform method of DataSet class to obtain the parameters of pixel coordinates to map coordinates, and then solve the binary quadratic equations to obtain pixel coordinates according to the formula of pixel coordinates to map coordinates and the specified map coordinate points. This calls the ReadRaster method.

3, Implemented code

1. Program code

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";
        //Register driver
        gdal.AllRegister();
        //Open file to get dataset
        Dataset dataset = gdal.Open(filepath,
                gdalconstConstants.GA_ReadOnly);
        if (dataset == null) {
            System.out.println("open"+filepath+"fail"+gdal.GetLastErrorMsg());
            System.exit(1);
        }
        //Get driver
        Driver driver = dataset.GetDriver();
        //Get driver information
        System.out.println("driver long name: " + driver.getLongName());
        //Gets the number of rasters
        int bandCount = dataset.getRasterCount();
        System.out.println("RasterCount: " + bandCount);
        //Construct an affine transformation parameter array and obtain data
        double[] gt = new double[6];
        dataset.GetGeoTransform(gt);
        System.out.println("Affine transformation parameters"+Arrays.toString(gt));

        //Specify latitude and longitude
        double Latitude  = 86.053;
        double longitude = 16.529;

        //Longitude and latitude are converted to grid pixel coordinates
        int[] ColRow=Coordinates2ColRow(gt,Latitude,longitude);

        //Determine whether the coordinates are out of range
        if(ColRow[0]<0||ColRow[1]<0||ColRow[0]>dataset.getRasterXSize()||ColRow[1]>dataset.getRasterYSize()){
            System.out.println(Arrays.toString(ColRow)+"Coordinate value is out of grid range!");
            return;
        }

        //Traverse the band, obtain the value of each band corresponding to the point and print it to the screen
        for (int i = 0;i<bandCount;i++){
            Band band = dataset.GetRasterBand(1);
            double[] values = new double[1];
            band.ReadRaster(ColRow[0], ColRow[1], 1, 1, values);
            double value = values[0];
            System.out.println("Abscissa:"+Latitude +","+"Ordinate:"+longitude);
            System.out.format("Band"+(i+1)+": %s", value);
        }

    }

    /**
     * Convert map coordinates to grid pixel coordinates
     * @param gt Affine transformation parameters
     * @param X Abscissa
     * @param Y Ordinate
     * @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];
//        Elimination method for solving binary linear equations
//        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])

        //Round down. If you round up, the calculation result will be too large, and the data of adjacent pixels will be read later
        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. Operation results

3. Verify

Here, QGIS is used to open the data file, and the recognition function is used to verify the values of relevant points.

  It can be found that the program analysis result is consistent with the QGIS analysis result.

4, Supplementary notes

The meaning of gt [] obtained by dataset.GetGeoTransform is as follows:

//gt[0]   x coordinate of upper left corner
//gt[1]   East West resolution
//gt[2]   Rotation angle 1, 0 indicates that the image is "north up"
//gt[3]   y coordinate of upper left corner
//gt[4]   Rotation angle 2, 0 indicates that the image is "north up"
//gt[5]   East West resolution

Map coordinate x=    Upper left corner x coordinate + number of columns * East-West resolution + number of rows * rotation angle 1;

Map coordinate y=    Y coordinate of upper left corner + number of columns * rotation angle 2  + Number of rows * resolution in north-south direction;

The input parameters of band.ReadRaster have the following meanings: 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 abscissa of partial origin position to be read

yoff wants to read the vertical image coordinates of the partial origin position

xsize specifies the side length of the rectangle to read part of the image.  

ysize specifies the rectangular side length of the part of the image to be read.  

Array is an array used to store the read data.

The data of this paper is made by the author, which has no business significance. If necessary, you can download: Demo.tif 

Topics: Java gis