一展鸡、背景介紹

由于項目中需要用到大范圍tiff的圖像（全中國30米分辨率的dem影像），并且需要單點獲取高程埃难，以及實現(xiàn)部分范圍的dem裁切與獲取趨于范圍極值莹弊，當時在網(wǎng)上查找的部分，很多都不滿足預期涡尘，或者計算結(jié)果與實際并不夠契合忍弛，因此單獨開一篇專門講這塊內(nèi)容。

二考抄、 依賴引入

引入基礎(chǔ)的依賴包细疚，使用中間組件完成所需功能。

       <!-- 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>

由于組件之間存在交互關(guān)系川梅，因此將其全部放到一起疯兼。

三、高程操作

3.0 pre-step: 基礎(chǔ)環(huán)境配置

由于接下去的操作都需要依賴同一份數(shù)據(jù)源挑势，因此在最開始的時候镇防，需要先封裝關(guān)于數(shù)據(jù)源的相關(guān)操作。

3.0.1 基礎(chǔ)數(shù)據(jù)源初始化

3.0.1.1 GridCoverage2D創(chuàng)建

public void gridCoverage2D() throws IOException {
        //讀取基礎(chǔ)dem文件
        File file = new File(demFilePath);
        //加載文件所需基礎(chǔ)參數(shù)與投影參數(shù)
        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潮饱，獲得基礎(chǔ)GridCoverage2D 来氧。
        GridCoverage2D coverage = tifReader.read(null);
}

3.0.1.2 TIFFImage創(chuàng)建

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

3.1 獲取指定經(jīng)緯度高程

通過獲取gridCoverage2D 的的投影信息，加載投影并使用經(jīng)緯度進行計算位置點高程

    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 獲取區(qū)域范圍極值

 @ApiImplicitParams({
            @ApiImplicitParam(name = "geometry", value = "飛行區(qū)域的數(shù)據(jù)", required = true, paramType = "GeoJSON"),
            @ApiImplicitParam(name = "geometry_type", value = "圖形數(shù)據(jù)格式(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
    /**
     * 獲取區(qū)域的高程點列表
     * @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);
    }


    /**
     * 計算數(shù)組極值
     *
     * @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 區(qū)域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給出一個基礎(chǔ)圖層，但這個基礎(chǔ)圖層使用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];
            //經(jīng)緯度轉(zhuǎn)點楔壤，X是經(jīng)度，表寬惯驼；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]);
                }
            }
            /* 更新參數(shù) */
            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();
        }
    }


    /**
     * 根據(jù)幾何模型進行影像切割
     *
     * @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);
    }

四、總結(jié)

總體來說祟牲，整個的實現(xiàn)流程基本上可以概括為[引入依賴]->[加載數(shù)據(jù)源]->[業(yè)務(wù)數(shù)據(jù)讀取]隙畜。讀取DEM的高程流程相對簡單，獲取范圍內(nèi)的高程極值操作也不復雜说贝，主要是需要算出范圍對應的格網(wǎng)议惰，以及獲取每個格網(wǎng)像元的高程值；較為復雜的是DEM的圖像導出乡恕，由于使用geotools導出的dem圖像存在一定的偏移與投影問題换淆，導致其他平臺軟件無法正常讀取dem的一些屬性信息，因此在操作上几颜，先使用geotools自定義的導出倍试，先輸出一份DEM，再使用正常平臺導出的dem文件蛋哭，通過讀取它的屬性信息县习，賦值給空的dem影像，再將像元的高程值寫入谆趾，最后輸出成一份符合規(guī)范的dem躁愿。此操作較為繁瑣，但基礎(chǔ)滿足功能所需沪蓬，此為拋磚引玉彤钟，如有更好的辦法實現(xiàn)，望諸位不吝賜教跷叉，將感激不盡逸雹！