GDAL，我暗恋你很久了！

2021-11-06  本文已影响0人  WongBynn

初识GDAL

GDAL（Geospatial Data Abstraction Library），是一个用于矢量和栅格格式数据之间的转换器。由开源地理空间基金会在 X/MIT 风格的开源许可下发布。作为一个库，它为所有支持的格式向调用应用程序提供单个栅格抽象数据模型和单个矢量抽象数据模型。它还带有各种有用的命令行实用程序，用于数据转换和处理。该新闻页面描述了2021年10月GDAL / OGR 3.3.3版本。

GIS的基础设施-GDAL

个人认为，GDAL可以称为GIS的基础设施 ，为什么这么说呢？因为GDAL无论是在开源GIS软件还是商业GIS软件中都得到了广泛的应用。我们通过下表来看看，哪些软件都在使用GDAL。

软件名称 用途
ESRI ArcGIS 9.2+ A popular GIS platform.
ERDAS ER Viewer Image viewer for very large JPEG 2000 and ECW files.
GeoServer a open source software server written in Java that allows users to share and edit geospatial data.
Google Earth A 3D world viewer.
IDRISI A GIS and Image Processing Windows Desktop application. Uses GDAL to import/export/warp raster data.
NASA Ames Stereo Pipeline Software for creating terrain models and ortho images from planetary stereo images.
QGIS A cross platform desktop GIS.
SkylineGlobe The Skyline suite of interactive applications allows you to build, view, query and analyze customized, virtual 3D landscapes.
PostgreSQL OGR Foreign Data Wrapper Expose OGR layer as PostgreSQL foreign tables.
PostGIS spatial database extender for PostgreSQL: The raster loader and many of the raster SQL functions rely on GDAL.
PDAL Point Cloud Data Abstraction Library
GRASS GIS A raster/vector open source GIS that uses GDAL for raster/vector import and export (via r.in.gdal/r.out.gdal)
NASA WorldWind Multiplatform virtual globe library to quickly and easily create interactive visualizations of 3D globes, map and geographical information.
R A free software environment for statistical computing and graphics, with bindings to GDAL via the rgdal package.
WindNinja wind model for fire behavior modeling.

知道大家英语好，就不翻译了哦！单词都很简单~

GDAL与OGR的关系

OGR曾经是一个独立的矢量IO库，灵感来自于OpenGIS Simple Features，它是从GDAL中分离出来的。在GDAL 2.0版本中，GDAL和OGR组件被集成在一起。

OGR过去是OpenGIS简单特性参考实现的缩写。然而，由于OGR不完全符合OpenGIS简单特性规范，并且没有被批准作为该规范的参考实现，因此名称改为OGR简单特性库。这个名字中OGR的唯一含义是历史性的。在库的源代码中，OGR也是用于类名、文件名等的前缀。

搭建GDAL开发环境

1. 搭建GDAL的Windows开发环境

由于在实际开发中，我们大多都是使用Windows机器来进行工作，因此搭建windows的开发环境就显得格外重要。但是在windows平台下编译GDAL步骤较为繁琐，且存在一定的坑，对于非C/C++/Qt开发人员来说需要额外的精力去学习和摸索。如果你感兴趣的话，可以参考这篇博主的文章完成GDAL的编译。

为了节省时间，并且能够多快好省的直接进入GDAL的学习环节，我们这里使用别人编译好的GDAL jar包和动态链接库文件。

点我进入下载编译好的GDAL页面

下载的页面会跳转至GISInternals页面，关于这个网站的介绍如下所示：

About GISInternals

GISInternals is an online system for creating daily built binary packages for the GDAL and MapServer projects. The anchestor of this system has been provided to be the Windows buildslaves for the osgeo buildbot in 2007. The build system in the current form (providing downloadable packages) has been set up in 2009. As of this time, the system has been continuously improved by adding more and more packages to make the life easier for the users and developers of these open source projects. During this time the number of the visitors and the amount of the downloads continues to grow, and the site has been visited from more than 160 countries around the world.

GISInternals 是一个在线系统，用于为 GDAL 和 MapServer 项目创建每日构建的二进制包。 2007年已经提供了这个系统的源头，作为osgeo buildbot的Windows buildslaves。目前形式的构建系统（提供可下载的包）已经在2009年建立起来了。截至目前，该系统一直在不断改进 通过添加越来越多的包，让这些开源项目的用户和开发人员的生活更轻松。 在此期间，访问者数量和下载量持续增长，该站点已被全球 160 多个国家/地区访问。

由此可见GISInternals网站是一个专门用于构建GDAL和MapServer工具网站。

进入下载页面选择合适的版本，我们这里选择release-1900-x64-gdal-2-3-1-mapserver-7-2-0.zip的版本，下载好之后，解压文件后文件列表如下所示：

image.png

至此，windows环境下的gdal 环境搭建就完成了，稍后我们会进行实战操作···

2. 搭建GDAL 的CentOS运行环境

我们编写好的应用程序，经过严格的测试之后都会上线运行在CentOS服务器上，当然了，你可能会使用其它类型的OS服务器，但这里我们一CentOS 7版本为例，因为centOS更有代表性。

你可以在公众号后台私信 gdal，获取在CentOS 7平台编译gdal需要的依赖包。依赖包如下图：

image.png
如果这个版本你觉得比较低，你也可以去这里下载你想过要的gdal源码包版本，然后去寻找对应版本的geos和proj依赖包。如果这里你有任何疑惑，可以私信我交流~

在你开始进行编译前得先确保服务器具有gcc-c++、gcc、swig、ant的 基础环境，如果没有那就使用yum 工具在线安装即可。如果是离线安装，依赖环境需要从安装镜像里面去下载安装。你需要把安装进行挂载至服务器上，然后进行安装这些环境，如果你有任何疑问可以私信与我交流。

之后编译 gdal 的 java 库：

[root@wb gdal-2.1.3]#cd swig/java

[root@wb java]# make

编译完成后，在 “……/gdal-2.1.3/swig/java”目录 下生成 gdal.jar 和

libgdalconstjni.so、libgdaljni.so、libogrjni.so、libosrjni.so 等 4 个 so 文件。

拷贝 4 个 so 文件到 java.library.path 路径中。

GDAL开发实战

你可以选择新建一个SpringBoot工程,也可以使用你已有的工程，前提是你得完成前面我们提到的环境搭建的工作，确保我们能在SpringBoot工程中可以引用gdal的功能接口；

为了方便起见，我建议你使用已有的SpringBoot工程进行开始我们的学习，我自己就是在一个专门进行新技术学习的SpringBoot工程里面进行学习。

实战1 读取影像的元信息

在下面的代码中，我们对影响做了元信息的提取，并对影像的元信息进行了简单的处理和组装，并返回这些结果。

package com.bin.utils.gis;

import com.alibaba.fastjson.JSONObject;
import com.geovis.bin.utils.DateUtil;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;
import java.io.File;

/**
 * @Author Wangb
 * @Date 2021/11/6 22:08.
 */
public class GDALUtil {

    /**
     * 读取影像元信息
     *
     * @param fileName
     * @return
     */
    public static JSONObject parseImage(String fileName) {
        JSONObject result = new JSONObject();
        try {
            //gdal使用前的所必须要有的注册语句
            gdal.AllRegister();
            long start = System.currentTimeMillis();
            // 读取影像数据
            Dataset dataset = gdal.Open(fileName, gdalconstConstants.GA_ReadOnly);
            if (dataset == null) {
                System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
                System.err.println(gdal.GetLastErrorMsg());
                return null;
            }

            Driver driver = dataset.GetDriver();
            //System.out.println("Driver: " + driver.getShortName() + "/" + driver.getLongName());
            // 读取影像信息
            String satellite_type = "CAR-" + fileName.substring(103, 105);
            resultJp2.put("卫星名称",satellite_type);
            result.put("参考坐标系标识", dataset.GetProjectionRef());
            result.put("坐标系标识", dataset.GetProjection());
            result.put("类型", driver.getShortName());
            result.put("列", dataset.getRasterXSize());
            result.put("行", dataset.getRasterYSize());
            result.put("band", dataset.getRasterCount());
            result.put("pixelsType", dataset.GetRasterBand(1).getDataType());// 像素类型

            double[] transform = dataset.GetGeoTransform();
            result.put("MinX", transform[0]);
            result.put("MaxX", (transform[0] + dataset.getRasterXSize() * transform[1]));
            result.put("MinY", (transform[3] - dataset.getRasterYSize() * transform[1]));
            result.put("MaxY", transform[3]);
            result.put("resolution", transform[1]); // 分辨率

            // 卫星四至坐标
            String geom = String.format("Polygon((%s %s,%s %s,%s %s,%s %s ,%s %s ))",
                    Double.parseDouble(result.getString("MaxX")),
                    Double.parseDouble(result.getString("MinY")),
                    Double.parseDouble(result.getString("MaxX")),
                    Double.parseDouble(result.getString("MaxY")),
                    Double.parseDouble(result.getString("MinX")),
                    Double.parseDouble(result.getString("MaxY")),
                    Double.parseDouble(result.getString("MinX")),
                    Double.parseDouble(result.getString("MinY")),
                    Double.parseDouble(result.getString("MaxX")),
                    Double.parseDouble(result.getString("MinY"))

            );

            // 获取文件拍摄时间
            File file = new File(fileName);
            String shotTime = file.getName().substring(36, 40) + "-" + file.getName().substring(34, 36) + "-" + file.getName().substring(32, 34);

            // 获取卫星名称
            String satelliteName = file.getName().substring(0, 3);

            // 获取拍摄中心点坐标
            String center = String.format("POINT(%s %s)",
                    (transform[0] + (transform[0] - dataset.getRasterYSize() * transform[1])) / 2,
                    (transform[3] + (transform[3] - dataset.getRasterYSize() * transform[1])) / 2
            );
            // 获取文件格式
            String format = fileName.substring(fileName.lastIndexOf(".") + 1);

            result.put("center", center);
            result.put("geom", geom);

            result.put("shotTime", DateUtil.parseDate(shotTime));
            result.put("satelliteName", "CAR-" + satelliteName.substring(1, 3));
            result.put("format", format);

           String projection = dataset.GetProjection();
                                            System.out.println(gdalconstConstants.GDT_Byte);
            dataset.delete();
            long end = System.currentTimeMillis();
            System.out.println(end - start);
            gdal.GDALDestroyDriverManager();
            return result;

        } catch (Exception e) {
            e.printStackTrace();
        }

        return null;

    }

    public static void main(String[] args) {
        parseImage("D:\\01_资料\\GIS\\GF1_PMS1_E121.3_N39.7_20170429_L1A0002339797/GF1_PMS1_E121.3_N39.7_20170429_L1A0002339797-MSS1.tiff");
    }
}

获取影像元信息的程序运行结果如下图所示：


image.png

实战2 影像的裁切

package com.geovis.bin.utils.gis;

import com.alibaba.fastjson.JSONObject;
import com.geovis.bin.utils.DateUtil;
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;
import org.gdal.gdalconst.gdalconstConstants;

import java.io.File;

/**
 * @Author Wangb
 * @Date 2021/6/11 19:08.
 */
public class GDALUtil {
    public static synchronized void  cutImage(String fileName) {
        gdal.AllRegister();
        Dataset rds = gdal.Open(fileName, gdalconst.GA_ReadOnly);
        //宽、高、波段数
        int b = rds.getRasterCount();
        //从波段中获取影像的数据类型，gdal中波段索引从1开始
        int dataType = rds.GetRasterBand(1).GetRasterDataType();

        //六参数信息
        double[] geoTransform = rds.GetGeoTransform();

        /**
         * 设置要裁剪的起始像元位置，以及各方向的像元数
         * 这里表示从像元(5000, 5000)开始，x方向和y方向各裁剪1000个像元
         * 最终结果就是一幅1000*1000的影像
         */
        int startX = 5000;
        int startY = 5000;
        int clipX = 1000;
        int clipY = 1000;

        //计算裁剪后的左上角坐标
        geoTransform[0] = geoTransform[0] + startX * geoTransform[1] + startY * geoTransform[2];
        geoTransform[3] = geoTransform[3] + startX * geoTransform[4] + startY * geoTransform[5];

        //创建结果图像
        Driver driver = gdal.GetDriverByName("GTIFF");
        String file = "D:\\01_资料\\GIS\\test.tiff";
        Dataset outputDs = driver.Create(file, clipX, clipY, b, dataType);
        outputDs.SetGeoTransform(geoTransform);
        outputDs.SetProjection(rds.GetProjection());

        try {
            //按band读取
            for (int i = 0; i < clipY; i++) {
                //按行读取
                for (int j = 1; j <= b; j++) {
                    Band orgBand = rds.GetRasterBand(j);
                    int[] cache = new int[clipX];
                    //从位置x开始，只读一行
                    orgBand.ReadRaster(startX, startY + i, clipX, 1, cache);
                    Band desBand = outputDs.GetRasterBand(j);
                    //从左上角开始，只写一行
                    desBand.WriteRaster(0, i, clipX, 1, cache);
                    desBand.FlushCache();
                }
            }
        } catch (Exception e) {
            e.printStackTrace();
        } finally {

            //释放资源
            rds.delete();
            outputDs.delete();
        }

    }

    public static void main(String[] args) {
        String fileName = "D:\\01_资料\\GIS\\GF1_PMS1_E121.3_N39.7_20170429_L1A0002339797/GF1_PMS1_E121.3_N39.7_20170429_L1A0002339797-MSS1.tiff";

//        parseImage(fileName);
        cutImage(fileName);
    }
}

影像裁切的运行截图如下，虽然成功的进行了影像的裁切，但是控制台会爆红，原因暂时还不清楚，有知道的小伙伴可以私信一下哈~

image.png
值得注意的是不论你是获取影像的元信息还是进行影像的裁切，都应该对方法进行并发控制，最简单的实现就是加上synchronized关键字，这样做虽然影响了并发时性能，但是却避免了OOM异常，因为每一次读取影像，都会把整个tif读进内存。

点击跳转至gdal API手册

后面我会继续进行GDAL其它常用功能的学习，也希望这篇文章可以帮助到你···

上一篇下一篇

猜你喜欢

热点阅读