利用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
來源：CSDN
原文：https://blog.csdn.net/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
來源：簡書