GIS_GeoTools栅格数据分析之图像变化检测

桃扇骨 2022-11-07 13:30 229阅读 0赞

*  [功能需求：给定同一区域不同时间的无人机影像数据，求出区域内影像变化部分，并矢量化成GeoJSON返回给前端。][GeoJSON]
 *  [1.将两幅图像进行相减与二值化操作][1.]
 *  [2.调用geoTools的PolygonExtractionProcess将图像相减操作结果进行矢量化][2._geoTools_PolygonExtractionProcess]
 *  [3.对矢量化后的多边形对象进行过滤，删除面积过小的细碎多边形][3.]
 *  [4.将最终结果以GeoJSON格式返回][4._GeoJSON]

功能需求：给定同一区域不同时间的无人机影像数据，求出区域内影像变化部分，并矢量化成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);
            }
        }

**转载自：**https://blog.csdn.net/weixin\_34226706/article/details/87275194

### 相关推荐 ###

*  [Geo\_tools—–Eclipse快速入门 显示shp地图][Geo_tools_Eclipse_ _shp]
 *  [GEoTools学习笔记—显示三维效果–geotools-renderer3d][GEoTools_geotools-renderer3d]
 *  [JTS Geometry之间的关系][JTS Geometry]
 *  [geotools获取shapefile 属性名乱码问题][geotools_shapefile]
 *  [geotools][]
 *  [【Neo4j07】CQL命令\_Sorting排序][Neo4j07_CQL_Sorting]

[GeoJSON]: https://www.giserdqy.com/secdev/geotools/27924/#%E5%8A%9F%E8%83%BD%E9%9C%80%E6%B1%82%EF%BC%9A%E7%BB%99%E5%AE%9A%E5%90%8C%E4%B8%80%E5%8C%BA%E5%9F%9F%E4%B8%8D%E5%90%8C%E6%97%B6%E9%97%B4%E7%9A%84%E6%97%A0%E4%BA%BA%E6%9C%BA%E5%BD%B1%E5%83%8F%E6%95%B0%E6%8D%AE%EF%BC%8C%E6%B1%82%E5%87%BA%E5%8C%BA%E5%9F%9F%E5%86%85%E5%BD%B1%E5%83%8F%E5%8F%98%E5%8C%96%E9%83%A8%E5%88%86%EF%BC%8C%E5%B9%B6%E7%9F%A2%E9%87%8F%E5%8C%96%E6%88%90GeoJSON%E8%BF%94%E5%9B%9E%E7%BB%99%E5%89%8D%E7%AB%AF%E3%80%82
[1.]: https://www.giserdqy.com/secdev/geotools/27924/#1%E5%B0%86%E4%B8%A4%E5%B9%85%E5%9B%BE%E5%83%8F%E8%BF%9B%E8%A1%8C%E7%9B%B8%E5%87%8F%E4%B8%8E%E4%BA%8C%E5%80%BC%E5%8C%96%E6%93%8D%E4%BD%9C
[2._geoTools_PolygonExtractionProcess]: https://www.giserdqy.com/secdev/geotools/27924/#2%E8%B0%83%E7%94%A8geoTools%E7%9A%84PolygonExtractionProcess%E5%B0%86%E5%9B%BE%E5%83%8F%E7%9B%B8%E5%87%8F%E6%93%8D%E4%BD%9C%E7%BB%93%E6%9E%9C%E8%BF%9B%E8%A1%8C%E7%9F%A2%E9%87%8F%E5%8C%96
[3.]: https://www.giserdqy.com/secdev/geotools/27924/#3%E5%AF%B9%E7%9F%A2%E9%87%8F%E5%8C%96%E5%90%8E%E7%9A%84%E5%A4%9A%E8%BE%B9%E5%BD%A2%E5%AF%B9%E8%B1%A1%E8%BF%9B%E8%A1%8C%E8%BF%87%E6%BB%A4%EF%BC%8C%E5%88%A0%E9%99%A4%E9%9D%A2%E7%A7%AF%E8%BF%87%E5%B0%8F%E7%9A%84%E7%BB%86%E7%A2%8E%E5%A4%9A%E8%BE%B9%E5%BD%A2
[4._GeoJSON]: https://www.giserdqy.com/secdev/geotools/27924/#4%E5%B0%86%E6%9C%80%E7%BB%88%E7%BB%93%E6%9E%9C%E4%BB%A5GeoJSON%E6%A0%BC%E5%BC%8F%E8%BF%94%E5%9B%9E
[Geo_tools_Eclipse_ _shp]: https://www.giserdqy.com/secdev/geotools/28064/
[GEoTools_geotools-renderer3d]: https://www.giserdqy.com/secdev/geotools/28061/
[JTS Geometry]: https://www.giserdqy.com/secdev/geotools/28056/
[geotools_shapefile]: https://www.giserdqy.com/secdev/geotools/28058/
[geotools]: https://www.giserdqy.com/secdev/geotools/28052/
[Neo4j07_CQL_Sorting]: https://www.giserdqy.com/secdev/geotools/28048/