DEM地形操作(geotools方式与nga方式)

一、背景介绍

由于项目中需要用到大范围tiff的图像（全中国30米分辨率的dem影像），并且需要单点获取高程，以及实现部分范围的dem裁切与获取趋于范围极值，当时在网上查找的部分，很多都不满足预期，或者计算结果与实际并不够契合，因此单独开一篇专门讲这块内容。

二、 依赖引入

引入基础的依赖包，使用中间组件完成所需功能。


       <!-- geotools依赖 -->
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-shapefile</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-main</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-epsg-hsql</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-geotiff</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-geojson</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-api</artifactId>
            <version>20.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-metadata</artifactId>
            <version>22.0</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-opengis</artifactId>
            <version>22.0</version>
            <scope>compile</scope>
        </dependency>

        <!-- gdal依赖 -->
        <dependency>
            <groupId>org.gdal</groupId>
            <artifactId>gdal</artifactId>
            <version>3.0.0</version>
        </dependency>
        <!-- nga依赖 -->
        <dependency>
            <groupId>mil.nga</groupId>
            <artifactId>tiff</artifactId>
            <version>2.0.2</version>
        </dependency>
        <!-- jts依赖 -->
        <dependency>
            <groupId>org.locationtech.jts</groupId>
            <artifactId>jts-core</artifactId>
            <version>1.16.1</version>
        </dependency>

由于组件之间存在交互关系，因此将其全部放到一起。

三、高程操作

3.0 pre-step: 基础环境配置

由于接下去的操作都需要依赖同一份数据源，因此在最开始的时候，需要先封装关于数据源的相关操作。

3.0.1 基础数据源初始化

3.0.1.1 GridCoverage2D创建

public void gridCoverage2D() throws IOException {
        //读取基础dem文件
        File file = new File(demFilePath);
        //加载文件所需基础参数与投影参数
        Hints tiffHints = new Hints();
        tiffHints.add(new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE));
        tiffHints.add(new Hints(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, DefaultGeographicCRS.WGS84));
        //tiff读取器
        GeoTiffReader tifReader = new GeoTiffReader(file, tiffHints);
        //读取tiff，获得基础GridCoverage2D 。
        GridCoverage2D coverage = tifReader.read(null);
}

3.0.1.2 TIFFImage创建

    public TIFFImage tiffImage() {
        File file = new File(demFileMin);
        try {
            return TiffReader.readTiff(file);
        } catch (IOException e) {
            e.printStackTrace();
            return null;
        }
    }

3.1 获取指定经纬度高程

通过获取gridCoverage2D 的的投影信息，加载投影并使用经纬度进行计算位置点高程

    public double getHeight(double lon, double lat) {
        CoordinateReferenceSystem crs = gridCoverage2D.getCoordinateReferenceSystem2D();

        DirectPosition position = new DirectPosition2D(crs, lon, lat);
        int[] results = (int[]) gridCoverage2D.evaluate(position);
        results = gridCoverage2D.evaluate(position, results);
        return results[0];
    }

3.2 获取区域范围极值

 @ApiImplicitParams({
            @ApiImplicitParam(name = "geometry", value = "飞行区域的数据", required = true, paramType = "GeoJSON"),
            @ApiImplicitParam(name = "geometry_type", value = "图形数据格式(wkt,geojson)，目前支持wkt与geojson格式，默认wkt", required = true, paramType = "GeoJSON"),
    })
    public Result<?> ConventionalTrajectoryPlanning(
            @RequestParam(name = "geometry", required = true) String geometryString,
            @RequestParam(name = "geometry_type", required = true, defaultValue = "wkt") String geometryType,
            HttpServletRequest req) {
        GeometryJSON geometryJSON = new GeometryJSON();
        Geometry geometry = null;
        try {
            if (geometryType.equals("geojson")) {
                geometry = geometryJSON.readPolygon(new ByteArrayInputStream(geometryString.getBytes()));
            } else {
                geometry = reader.read(geometryString);
            }
            Geometry geometryMercator = JTS.transform(geometry, mathTransform);
            double area = geometryMercator.getArea() / 1000000;
            if (area > heightArea) {
                return Result.error("当前查询范围：" + String.format("%.2f", area) + "平方公里, " + "范围超出" + heightArea + "平方公里限制!");
            }

        } catch (IOException | ParseException | TransformException e) {
            e.printStackTrace();
            Result.error("GeoJSON输出异常，请检查输入的图形");
        }
        if (oConvertUtils.isEmpty(geometry)) return Result.error("输入图形为空，请检查");

        Geometry gridPolygon = getRegionGrid(geometry, TIFF_RESOLUTION);

        Coordinate[] coordinates = gridPolygon.getCoordinates();
        int arrLength = coordinates.length;
        double[] heightValueArr = new double[arrLength];
        for (int i = 0; i < arrLength; i++) {
            Coordinate coordinate = coordinates[i];
            heightValueArr[i] = geotoolsTerrainHeight.getHeight(coordinate.x, coordinate.y);
        }

        return Result.ok(maxAndMin(heightValueArr, arrLength));
    }

//引用的方法
//tiff的像元大小（30m*30m的DEM分辨率，一般可通过软件直接读取）
final double TIFF_RESOLUTION = 0.0002777777799991554275
    /**
     * 获取区域的高程点列表
     * @param geometry 输入的Geometry
     * @param tiffResolution tiff的原始分辨率
     * @return
     */
    private Geometry getRegionGrid(Geometry geometry, double tiffResolution) {

        GeometryFactory geometryFactory = new GeometryFactory()
        Envelope rect = geometry.getEnvelopeInternal();
        double height = rect.getHeight();
        double width = rect.getWidth();
        int numX = (int) Math.ceil(width / tiffResolution);
        int numY = (int) Math.ceil(height / tiffResolution);
        double dx = (width - numX * tiffResolution) / 2.0;
        double dy = (height - numY * tiffResolution) / 2.0;
        Geometry[][] nodes = new Geometry[numX][numY];
        for (int i = 0; i < numX; ++i) {
            for (int j = 0; j < numY; ++j) {
                double minX = dx + rect.getMinX() + i * tiffResolution;
                double minY = dy + rect.getMinY() + j * tiffResolution;
                double maxX = minX + tiffResolution;
                double maxY = minY + tiffResolution;
                Coordinate coord0 = new Coordinate(minX, minY);
                Coordinate coord2 = new Coordinate(maxX, minY);
                Coordinate coord3 = new Coordinate(maxX, maxY);
                Coordinate coord4 = new Coordinate(minX, maxY);
                Geometry box = geometryFactory.createPolygon(new Coordinate[]{coord0, coord2, coord3, coord4, coord0});
                if (box.intersects(geometry)) {
                    Geometry region = null;
                    try {
                        region = geometry.intersection(box);
                    } catch (Exception e) {
                        e.printStackTrace();
                        try {
                            region = box.intersection(geometry);
                        } catch (Exception ee) {
                            e.printStackTrace();
                            log.info("获取交点失败！box: " + box + "\r\n geometry:" + geometry);
                        }
                    }
                    nodes[i][j] = region;
                }
            }
        }
        List<Geometry> list = new ArrayList<Geometry>();
        for (int l = 0; l < numX; ++l) {
            for (int m = 0; m < numY; ++m) {
                Geometry region2 = nodes[l][m];
                if (region2 != null) {
                    list.add(region2);
                }
            }
        }
        return geometryFactory.buildGeometry(list);
    }


    /**
     * 计算数组极值
     *
     * @param a
     * @param length
     * @return
     */
    public JSONObject maxAndMin(double[] a, int length) {
        JSONObject jsonObject = new JSONObject();
        double Max, Min;
        double[] b, c;
        if (length % 2 == 0) {
            b = new double[length / 2];
            c = new double[length / 2];
        } else {
            b = new double[length / 2 + 1];
            c = new double[length / 2 + 1];
            b[length / 2] = a[length - 1];
            c[length / 2] = a[length - 1];
        }
        for (int i = 0, j = 0; i < length - 1; i += 2, j++) {
            if (a[i] >= a[i + 1]) {
                b[j] = a[i];
                c[j] = a[i + 1];
            } else {
                c[j] = a[i];
                b[j] = a[i + 1];
            }
        }

        Max = b[0];
        Min = c[0];
        for (int i = 1; i < b.length; i++) {
            if (Max < b[i]) Max = b[i];
        }
        for (int i = 1; i < c.length; i++) {
            if (Min > c[i]) Min = c[i];
        }
        jsonObject.put("max", Max);
        jsonObject.put("min", Min);
        jsonObject.put("avg", Double.valueOf(String.format("%.2f",Arrays.stream(a).average().orElse(Double.NaN))));

        return jsonObject;
    }

3.3 区域DEM文件导出


    @ApiOperation(value = "高程辅助-返回DEM文件流", notes = "高程辅助,返回DEM文件流")
    @ApiImplicitParams({
            @ApiImplicitParam(name = "wkt", value = "wkt格式的地理图形", required = true, paramType = "String"),
    })
    public void getDEM(
            @RequestParam(name = "wkt", required = true) String wkt,
            HttpServletRequest req, HttpServletResponse response) {

        try {
            /* 先用geotools给出一个基础图层，但这个基础图层使用nga的tiff包读不到像元尺度，因此再使用nga的tiff包重写 */
            Geometry geometry = reader.read(wkt);
            String prefix = UUID.randomUUID().toString().replace("-","");
           //定义临时文件的后缀
            String suffix =".tif";
            File tempOutFile = File.createTempFile(prefix,suffix);
            GridCoverage2D gridCoverage2D =geotoolsCoverageCrop(StringOperationUtils.gridCoverage2D, geometry);

            GeoTiffWriter writer = new GeoTiffWriter(tempOutFile, tiffHints);
            GeoTiffFormat geoTiffFormat = new GeoTiffFormat();
            ParameterValueGroup writeParameters = geoTiffFormat.getWriteParameters();
            java.util.List<GeneralParameterValue> valueList = writeParameters.values();
            writer.write(gridCoverage2D, valueList.toArray(new GeneralParameterValue[valueList.size()]));
            writer.dispose();

            /* 使用nga的tiff包重写 */
            TIFFImage tmpImage= TiffReader.readTiff(tempOutFile);
            FileDirectory tmpFileDirectory = tmpImage.getFileDirectory();
            Rasters rasters = tmpFileDirectory.readRasters();
            double[] lowerCorner = gridCoverage2D.getEnvelope().getLowerCorner().getCoordinate();
            double[] upperCorner = gridCoverage2D.getEnvelope().getUpperCorner().getCoordinate();
            //左上角定点
            double minX = lowerCorner[0];
            double maxY = upperCorner[1];
            //经纬度转点，X是经度，表宽；y是纬度，表高
            FileDirectory fileDirectory = tiffImage.getFileDirectory();
            //对新影像赋值
            int width = rasters.getWidth(), height = rasters.getHeight();
            Rasters newRaster = new Rasters(width, height, fileDirectory.getSamplesPerPixel(),
                    fileDirectory.getBitsPerSample().get(0), TiffConstants.SAMPLE_FORMAT_SIGNED_INT);
            for (int y = 0; y < height; y++) {
                for (int x = 0; x < width; x++) {
                    newRaster.setFirstPixelSample(x , y , rasters.getPixel(x, y)[0]);
                }
            }
            /* 更新参数 */
            int rowsPerStrip = newRaster.calculateRowsPerStrip(TiffConstants.PLANAR_CONFIGURATION_PLANAR);
            fileDirectory.setRowsPerStrip(rowsPerStrip);
            fileDirectory.setImageHeight(height);
            fileDirectory.setImageWidth(width);
            fileDirectory.setWriteRasters(newRaster);
            List<Double> doubles = new ArrayList<>();
            doubles.add(0.00);
            doubles.add(0.00);
            doubles.add(0.00);
            doubles.add(minX);
            doubles.add(maxY);
            doubles.add(0.00);
            FileDirectoryEntry fileDirectoryEntry = new FileDirectoryEntry(FieldTagType.ModelTiepoint, FieldType.DOUBLE, 6, doubles);
            fileDirectory.addEntry(fileDirectoryEntry);
            TIFFImage tiffImage = new TIFFImage();
            tiffImage.add(fileDirectory);
            byte[] tiffBytes = TiffWriter.writeTiffToBytes(tiffImage);
            response.getOutputStream().write(tiffBytes, 0, tiffBytes.length);
            response.getOutputStream().flush();
            /* 删除临时文件 */
            tempOutFile.delete();
        } catch (ParseException | IOException e) {
            e.printStackTrace();
        }
    }


    /**
     * 根据几何模型进行影像切割
     *
     * @param coverage2D 原始影像
     * @param geom       几何模型
     */
    public static GridCoverage2D geotoolsCoverageCrop(GridCoverage2D coverage2D, Geometry geom) {
        org.opengis.geometry.Envelope envelope = new Envelope2D();
        ((Envelope2D) envelope).setCoordinateReferenceSystem(DefaultGeographicCRS.WGS84);
        org.locationtech.jts.geom.Envelope jtsEnv = geom.getEnvelopeInternal();
        ((Envelope2D) envelope).height = jtsEnv.getHeight();
        ((Envelope2D) envelope).width = jtsEnv.getWidth();
        ((Envelope2D) envelope).x = jtsEnv.getMinX();
        ((Envelope2D) envelope).y = jtsEnv.getMinY();
        CoverageProcessor processor = CoverageProcessor.getInstance();

        final ParameterValueGroup param = processor.getOperation("CoverageCrop").getParameters();
        param.parameter("Source").setValue(coverage2D);
        param.parameter("Envelope").setValue(envelope);

        return (GridCoverage2D) processor.doOperation(param);
    }

四、总结

总体来说，整个的实现流程基本上可以概括为[引入依赖]->[加载数据源]->[业务数据读取]。读取DEM的高程流程相对简单，获取范围内的高程极值操作也不复杂，主要是需要算出范围对应的格网，以及获取每个格网像元的高程值；较为复杂的是DEM的图像导出，由于使用geotools导出的dem图像存在一定的偏移与投影问题，导致其他平台软件无法正常读取dem的一些属性信息，因此在操作上，先使用geotools自定义的导出，先输出一份DEM，再使用正常平台导出的dem文件，通过读取它的属性信息，赋值给空的dem影像，再将像元的高程值写入，最后输出成一份符合规范的dem。此操作较为繁琐，但基础满足功能所需，此为抛砖引玉，如有更好的办法实现，望诸位不吝赐教，将感激不尽！

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 194,911评论 5赞 460
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 82,014评论 2赞 371
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 142,129评论 0赞 320
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,283评论 1赞 264
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 61,159评论 4赞 357
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,161评论 1赞 272
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,565评论 3赞 382
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,251评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,531评论 1赞 292
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,619评论 2赞 310
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,383评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,255评论 3赞 313
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,624评论 3赞 299
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,916评论 0赞 17
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,199评论 1赞 250
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,553评论 2赞 342
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,756评论 2赞 335

推荐阅读更多精彩内容

  • 2019－11－28 周四 阴
    今天感恩节哎，感谢一直在我身边的亲朋好友。感恩相遇！感恩不离不弃。 中午开了第一次的党会，身份的转变要...
    迷月闪星情阅读 10,534评论 0赞 11
  • 2019.11.28下午市电视台
    彩排完，天已黑
    刘凯书法阅读 4,167评论 1赞 3
  • 表情管理
    表情是什么，我认为表情就是表现出来的情绪。表情可以传达很多信息。高兴了当然就笑了，难过就哭了。两者是相互影响密不可...
    Persistenc_6aea阅读 123,198评论 2赞 7