在Java GDAL环境中读取Tiff文件中的DEM数据

1. GDAL简介

GDAL全称是Geospatial Data Abstraction Library（地理空间数据抽象库），是使用C/C++语言编写的用于读写空间数据（包括栅格数据和矢量数据）的一套跨平台开源库，它利用抽象数据模型来表达所支持的各种文件格式，还使用一系列命令行来进行数据转换和处理。现有的大部分GIS或者遥感平台，不论是商业软件ArcGIS，ENVI还是开源软件GRASS，QGIS，都使用了GDAL作为底层构建库。

GDAL最初是由Frank Warmerdam于1998年开始开发的，在GDAL1.3.2版本之后，正式有开源空间信息基金会(Open Source Geospatial Foundation，OSGeo)下的GDAL/OGR项目委员会对其进行维护。

GDAL库由OGR和GDAL项目合并而来，OGR主要用于空间要素矢量矢量数据的解析，GDAL主要用于空间栅格数据的读写。此外，空间参考及其投影转换使用开源库 PROJ.4进行。

目前，GDAL主要提供了三大类数据的支持：栅格数据，矢量数据以及空间网络数据（Geographic Network Model）。栅格数据格式包括Arc/Info ASCII Grid（ASC），GeoTiff（Tiff），Erdas Imagine Images（Img），ASCII DEM（DEM）等。矢量数据格式主要包括ESRI Shapefiles、S-57、SDTS、PostGIS、Oracle Spatial、Mapinfo mid/mif和Mapinfo TAB等。

GDAL提供了C/C++借口，并且通过SWIG提供了Python，Java，C#等的调用借口。当我们在Java中调用GDAL的API函数时，其实底层执行的是C/C++编译的二进制文件。

2. 由Tiff文件中读取DEM数据

下面以GDAL演示工程为例，介绍一下由Tiff文件中读取DEM数据的流程

1）使用GDAL的第一步，需要注册文件格式

import org.gdal.gdal.gdal;
gdal.AllRegister();

或者

org.gdal.gdal.gdal.AllRegister();

2）定义数据库，并以只读方式打开数据文件

Dataset poDataset = gdal.Open("srtm_60_05.tif", gdalconst.GA_ReadOnly);

"srtm_60_05.tif"是文件名，该文件名可以是工程目录下的文件，也可以是全路径的文件名。

此时，Tiff文件中的数据已经读到了poDataset中，后续的数据操作将围绕poDataset进行。

3）定义驱动接口

Driver hDriver = poDataset.GetDriver();

4）定义相应的数据波段

SRTM3的Tiff文件中只有一个波段，定义相应的数据波代码为：

org.gdal.gdal.Band poBand = poDataset.GetRasterBand(1);

若提取数据中的波段数目，可以用命令：poDataset.getRasterCount()

5）读取相应的数据

poBand.ReadRaster(x0, y0, x_points, y_points, values);

x0与y0为预读取数据在数据栅格中左上的行号与列号，x_points与y_points分别为沿x向与y向的点数，而values中将存放相应读出的数组，注意values应该是一维数组。

3. 其它相关命令

1） 数据驱动接口

Driver hDriver = poDataset.GetDriver();

2）文件格式

hDriver.GetDescription()

3）数据栅格的尺度

poDataset.getRasterXSize()；poDataset.getRasterYSize()

4）输出图像的投影信息

poDataset.GetProjectionRef()

如：GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]

5）获取仿射变换参数

double[] gt = new double[6];

poDataset.GetGeoTransform(gt);

在“上北下南”的图像中，系数gt(2)与gt(4) 为0， 系数gt(1) 是像素宽度，gt(5)是像素高度，(gt(0),gt(3))位置是栅格左上角像素的经纬度。

栅格坐标和地理坐标之间的仿射变换关系为

表示实际的地址坐标，即经纬度；为栅格数据中的坐标。

6）获取波段数据尺度

poBand.getXSize()；poBand.getYSize()

4. 示例

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

public class ReadDEM {

	public static void main(String[] args) {
		// 注册文件格式
		gdal.AllRegister();
		// 使用只读的方式打开图像
		Dataset poDataset = gdal.Open("srtm_60_05.tif", gdalconst.GA_ReadOnly);
		
		Driver hDriver = poDataset.GetDriver();
		
		// 定义相应的数据波段
		Band poBand = poDataset.GetRasterBand(1);
		
		// 读取相应的数据
		int values[] = new int[5];
		
		for (int x = 0; x < 4; x ++) {
			poBand.ReadRaster(10 + x, 20, 1, 5, values);
			for (int y = 0; y < 5; y++) {
				System.out.print(values[y] + ", ");
			}
			System.out.println("\n");
		}

		poDataset.delete();
			
	}

}

代码运行结果：

 通过调整ReadRaster的参数，就可以实现对栅格数据集中的特定数据进行读取。