利用Geotools来转换影像的坐标系[转]

使用Operations类下的resample方法(重采样)可以解决这个问题，它的方法定义如下： 

Coverage org.geotools.coverage.processing.Operations.resample(Coverage source, CoordinateReferenceSystem crs) throws CoverageProcessingException

 

File file = new File("xxxx.tif");
if(file.exists()){
Reader br = new Reader();
GridCoverage2D old2D = br.getGridCoverage2D(file);
final CoordinateReferenceSystem WGS = CRS.decode("EPSG:3857");
final CoordinateReferenceSystem sourceCRS = old2D.getCoordinateReferenceSystem();
System.out.println(String.format("源坐标系为: %s", sourceCRS.getName()));
GridCoverage2D new2D = (GridCoverage2D) Operations.DEFAULT.resample(old2D, WGS);
System.err.println(String.format("目标坐标系为: %s", new2D.getCoordinateReferenceSystem().getName()))；
}
---------------------
作者：WilsonOnIsland
来源：优快云
原文：https://blog.youkuaiyun.com/u013323965/article/details/79453383

 

1.计算指定像素坐标点的坡度

public float calcSlope(int cellX, int cellY, PlanarImage image) throws IOException {
        DecimalFormat df = new DecimalFormat("#.0000");
        final int[] dest = null;
        int e = image.getTile(image.XToTileX(cellX), image.YToTileY(cellY)).getPixel(cellX, cellY, dest)[0];
        int e1 = image.getTile(image.XToTileX(cellX - 1), image.YToTileY(cellY)).getPixel(cellX - 1, cellY, dest)[0];
        int e2 = image.getTile(image.XToTileX(cellX), image.YToTileY(cellY - 1)).getPixel(cellX, cellY - 1, dest)[0];
        int e3 = image.getTile(image.XToTileX(cellX + 1), image.YToTileY(cellY)).getPixel(cellX + 1, cellY, dest)[0];
        int e4 = image.getTile(image.XToTileX(cellX), image.YToTileY(cellY + 1)).getPixel(cellX, cellY + 1, dest)[0];
        int e5 = image.getTile(image.XToTileX(cellX - 1), image.YToTileY(cellY - 1)).getPixel(cellX - 1, cellY - 1,
                dest)[0];
        int e6 = image.getTile(image.XToTileX(cellX + 1), image.YToTileY(cellY - 1)).getPixel(cellX + 1, cellY - 1,
                dest)[0];
        int e7 = image.getTile(image.XToTileX(cellX + 1), image.YToTileY(cellY + 1)).getPixel(cellX + 1, cellY + 1,
                dest)[0];
        int e8 = image.getTile(image.XToTileX(cellX - 1), image.YToTileY(cellY + 1)).getPixel(cellX - 1, cellY + 1,
                dest)[0];
        double slopeWE = ((e8 + 2 * e1 + e5) - (e7 + 2 * e3 + e6)) / (8 * 2041.823085);// 东西方向坡度
        double slopeNW = ((e7 + 2 * e4 + e8) - (e6 + 2 * e2 + e5)) / (8 * 2041.823085);// 南北方向坡度
        double slope = 100*(Math.sqrt(Math.pow(slopeWE, 2) + Math.pow(slopeNW, 2)));
        return Float.parseFloat(df.format(slope));
    }

2.DEM数据测试用例

 @Test
    public void testCalcSlope() throws NoSuchAuthorityCodeException, FactoryException, IOException {
        String path = "D:\\workData\\geotiff\\testTiff.tif";
        String outputPath = "D:\\workData\\geotiff\\output.tif";
        File file = new File(path);
        // 设置tiff影像默认设置
        Hints tiffHints = new Hints();
        tiffHints.add(new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE));
        // 默认坐标系EPSG:3857 
//      tiffHints.add(new Hints(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, CRS.decode("EPSG:4326")));
        tiffHints.add(new Hints(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, DefaultGeographicCRS.WGS84));

        GeoTiffReader reader = new GeoTiffReader(file, tiffHints);

        GridCoverage2D coverage = reader.read(null);
        Envelope env = coverage.getEnvelope();
        PlanarImage image = (PlanarImage) coverage.getRenderedImage();
        int width = image.getWidth(); // Image Width
        int height = image.getHeight(); // Image Height

        // 计算每个栅格的坡度
        float[][] slopeData = new float[height][width];
        for (int i = 1; i < height + 1; i++) {
            for (int j = 1; j < width + 1; j++) {
                float slope = SlopeUtil.INSTANCE.calcSlope(j, i, image);
                slopeData[i - 1][j - 1] = slope;
            }
        }

        GridCoverageFactory factory = new GridCoverageFactory();
        GridCoverage2D outputCoverage = factory.create("test", slopeData, env);
        GeoTiffWriter writer = new GeoTiffWriter(new File(outputPath));
        writer.write(outputCoverage, null);
        writer.dispose();
    }

功能需求：给定同一区域不同时间的无人机影像数据，求出区域内影像变化部分，并矢量化成GeoJSON返回给前端。

1.将两幅图像进行相减与二值化操作

    public GridCoverage2D tiffSubtract(String sourceTiffPath, String targetTiffPath, float diffLimit)
            throws IOException {
        File sourceTiff = new File(sourceTiffPath);
        File targetTiff = new File(targetTiffPath);

        if (!sourceTiff.exists() || !targetTiff.exists()) {
            throw new FileNotFoundException(sourceTiffPath + " or " + targetTiffPath + " do not exist!");
        }

        // 中间数据tiff存储路径
        String tempTiff = sourceTiff.getParent() + File.separator + "output.tiff";
        // tiff文件坐标系设置
        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));

        GeoTiffReader sourceReader = new GeoTiffReader(sourceTiff, tiffHints);
        GeoTiffReader targetReader = new GeoTiffReader(targetTiff, tiffHints);

        GridCoverage2D sourceCoverage = sourceReader.read(null);
        GridCoverage2D targetCoverage = targetReader.read(null);
        RenderedImage sourceImage = sourceCoverage.getRenderableImage(0, 1).createDefaultRendering();
        RenderedImage targetImage = targetCoverage.getRenderableImage(0, 1).createDefaultRendering();
        Raster sourceRaster = sourceImage.getData();
        Raster targetRaster = targetImage.getData();

        int width = sourceRaster.getWidth();
        int height = sourceRaster.getHeight();
        // System.out.println("pixels : width:" + width + ";height:" + height);

        Envelope2D sourceEnv = sourceCoverage.getEnvelope2D();

        float[][] difference = new float[height][width];
        float s;
        float t;
        // 修改算法,提取差异值大于阈值的部分
        // 将图像二值化
        for (int x = 0; x < width - 1; x++) {
            for (int y = 0; y < height - 1; y++) {
//              System.out.println("x:" + x + ";y:" + y);
                s = sourceRaster.getSampleFloat(x, y, 1);
                t = targetRaster.getSampleFloat(x, y, 1);
                float diff = t - s;
                if (diff > diffLimit) {
                    difference[y][x] = 100f;
                } else {
                    difference[y][x] = 0f;
                }
            }
        }

        GridCoverageFactory factory = new GridCoverageFactory();
        GridCoverage2D outputCoverage = factory.create("subtractTiff", difference, sourceEnv);
        GeoTiffWriter writer = new GeoTiffWriter(new File(tempTiff));
        writer.write(outputCoverage, null);
        writer.dispose();
        return outputCoverage;
    }

2.调用geoTools的PolygonExtractionProcess将图像相减操作结果进行矢量化

    public String polygonExtraction(GridCoverage2D tiffCoverage, String shpPath)
            throws MismatchedDimensionException, IndexOutOfBoundsException, NoSuchAuthorityCodeException,
            ParseException, FactoryException, TransformException, SchemaException, IOException {
        PolygonExtractionProcess process = new PolygonExtractionProcess();
        SimpleFeatureCollection features = process.execute(tiffCoverage, 0, Boolean.TRUE, null, null, null, null);

        features = this.polygonPostprocess(features, 10d);

        SimpleFeatureType type = features.getSchema();

//      ShapeFileWriter.INSTANCE.write(shpPath, features, type);
        this.toGeoJSON(features);

        return shpPath;
    }

3.对矢量化后的多边形对象进行过滤，删除面积过小的细碎多边形

private SimpleFeatureCollection polygonPostprocess(SimpleFeatureCollection features, double aeraLimit)
            throws IndexOutOfBoundsException, ParseException, NoSuchAuthorityCodeException, FactoryException,
            MismatchedDimensionException, TransformException, SchemaException {
        //坐标转换,从4326转成3857
        CoordinateReferenceSystem dataCRS = DefaultGeographicCRS.WGS84;
        CoordinateReferenceSystem targerCRS = CRS.decode("EPSG:3857");
        boolean lenient = true; // allow for some error due to different datums
        MathTransform transform = CRS.findMathTransform(dataCRS, targerCRS, lenient);

        final SimpleFeatureType TYPE = DataUtilities.createType("Location",
                "the_geom:Polygon:srid=3857,DN:String,Aera:Double");

        List<SimpleFeature> projectedFeatureList = new ArrayList<SimpleFeature>();

        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
        WKTReader reader = new WKTReader(geometryFactory);
        SimpleFeatureBuilder builder = new SimpleFeatureBuilder(TYPE);

        SimpleFeatureIterator iterator = features.features();
        try {
            while (iterator.hasNext()) {
                SimpleFeature feature = iterator.next();
                Polygon polygon = (Polygon) feature.getDefaultGeometry();
                polygon = (Polygon) JTS.transform(polygon, transform);
                double aera = polygon.getArea();
                // 多边形面积大于阈值
                if (aera >= aeraLimit) {
                    builder.add(polygon);
                    builder.add(feature.getAttribute(1).toString());
                    builder.add(aera);
                    SimpleFeature tempFeature = builder.buildFeature(null);
                    projectedFeatureList.add(tempFeature);
                }
            }
        } finally {
            iterator.close();
        }

        System.out.println(projectedFeatureList.size());
        return new ListFeatureCollection(TYPE, projectedFeatureList);
    }

4.将最终结果以GeoJSON格式返回

    private void toGeoJSON(SimpleFeatureCollection featureCollection) {
        SimpleFeatureIterator it = featureCollection.features();
        GeoJsonWriter geoJsonWriter = new GeoJsonWriter();
        
        while(it.hasNext()) {
            SimpleFeature tempFeature = it.next();
            
            Geometry geometry = (Geometry) tempFeature.getDefaultGeometry();
            
            String json = geoJsonWriter.write(geometry);
            
            System.out.println(json);
        }
    }

 

作者：z362831561
链接：https://www.jianshu.com/p/6b7addbe8739
来源：简书
 

转载于:https://www.cnblogs.com/zany-hui/articles/11105923.html