Wechat: yu389741| Email: gisdqy@163.com

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

##### 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();

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();
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) {
SimpleFeature tempFeature = builder.buildFeature(null);
}
}
} 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);
}
}
``````