Coordinate Reference System Transformation

To change the coordinate reference system (CRS) of a raster image, the resample method from the Operations class is utilized. This method performs a reprojection of the source coverage to the target CRS.

The following example demonstrates how to reproject a GeoTIFF file to EPSG:3857 (Web Mercator):

public GridCoverage2D reprojectRaster(File inputFile) throws Exception {
    // Initialize reader with hints to ensure longitude is the first axis
    Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
    GeoTiffReader reader = new GeoTiffReader(inputFile, hints);
    
    GridCoverage2D sourceCoverage = reader.read(null);
    CoordinateReferenceSystem sourceCRS = sourceCoverage.getCoordinateReferenceSystem();
    System.out.println("Source CRS: " + sourceCRS.getName());

    // Define target CRS
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3857");

    // Perform the resampling operation
    GridCoverage2D reprojectedCoverage = (GridCoverage2D) Operations.DEFAULT.resample(sourceCoverage, targetCRS);
    System.out.println("Target CRS: " + reprojectedCoverage.getCoordinateReferenceSystem().getName());
    
    return reprojectedCoverage;
}

Calculating Slope from DEM Data

Slope calculation involves analyzing the elevation differences between a central pixel and its surrounding neighbors. The algorithm typically uses a 3x3 window to compute the rate of change in the x (West-East) and y (North-South) directions.

public float calculateSlope(int col, int row, PlanarImage image, double cellSize) {
    // Retrieve elevation values for the 3x3 neighborhood
    double northWest = getPixelElevation(image, col - 1, row - 1);
    double north     = getPixelElevation(image, col,     row - 1);
    double northEast = getPixelElevation(image, col + 1, row - 1);
    double west      = getPixelElevation(image, col - 1, row);
    double east      = getPixelElevation(image, col + 1, row);
    double southWest = getPixelElevation(image, col - 1, row + 1);
    double south     = getPixelElevation(image, col,     row + 1);
    double southEast = getPixelElevation(image, col + 1, row + 1);

    // Calculate gradients using Horn's formula
    double dx = ((northWest + 2 * west + southWest) - (northEast + 2 * east + southEast)) / (8 * cellSize);
    double dy = ((northWest + 2 * north + northEast) - (southWest + 2 * south + southEast)) / (8 * cellSize);

    // Calculate slope percentage
    double slopeRatio = Math.sqrt(dx * dx + dy * dy);
    return (float) (slopeRatio * 100.0);
}

private int getPixelElevation(PlanarImage image, int x, int y) {
    // Safely retrieve pixel value, handling boundary checks if necessary
    return image.getTile(image.XToTileX(x), image.YToTileY(y)).getSample(x, y, 0);
}

Change Detection and Vectorization

This process identifies changes between two temporal raster images (e.g., drone imagery) and converts the results into vector polygons. The workflow involves image subtraction, thresholding, polygon extraction, and filtering.

1. Image Subtraction and Binarization

First, the pixel values of the two images are compared. A threshold is applied to create a binary mask where significant differences are highlighted.

public GridCoverage2D createDifferenceMask(File sourceFile, File targetFile, float threshold) throws IOException {
    GeoTiffReader sourceReader = new GeoTiffReader(sourceFile);
    GeoTiffReader targetReader = new GeoTiffReader(targetFile);
    
    GridCoverage2D sourceCoverage = sourceReader.read(null);
    GridCoverage2D targetCoverage = targetReader.read(null);
    
    RenderedImage sourceImage = sourceCoverage.getRenderedImage();
    RenderedImage targetImage = targetCoverage.getRenderedImage();
    
    int width = sourceImage.getWidth();
    int height = sourceImage.getHeight();
    float[][] differenceData = new float[height][width];
    
    Raster sourceRaster = sourceImage.getData();
    Raster targetRaster = targetImage.getData();
    
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            float val1 = sourceRaster.getSampleFloat(x, y, 0);
            float val2 = targetRaster.getSampleFloat(x, y, 0);
            
            // Apply threshold to determine change
            if (Math.abs(val2 - val1) > threshold) {
                differenceData[y][x] = 100f; // Mark as changed
            } else {
                differenceData[y][x] = 0f;   // No change
            }
        }
    }
    
    GridCoverageFactory factory = new GridCoverageFactory();
    return factory.create("difference", differenceData, sourceCoverage.getEnvelope());
}

2. Polygon Extraction and Filtering

The binary raster is converted into vector polygons using PolygonExtractionProcess. Post-processing removes small, insignificant polygons based on an area threshold.

public void extractAndExportChangePolygons(GridCoverage2D diffCoverage, double minArea) throws Exception {
    // Convert raster to vector polygons
    PolygonExtractionProcess extractor = new PolygonExtractionProcess();
    SimpleFeatureCollection rawFeatures = extractor.execute(diffCoverage, 0, true, null, null, null, null);
    
    // Reproject to a metric CRS (e.g., EPSG:3857) for accurate area calculation
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:3857");
    MathTransform transform = CRS.findMathTransform(DefaultGeographicCRS.WGS84, targetCRS, true);
    
    SimpleFeatureTypeBuilder typeBuilder = new SimpleFeatureTypeBuilder();
    typeBuilder.setName("ChangePolygons");
    typeBuilder.setCRS(targetCRS);
    typeBuilder.add("the_geom", Polygon.class);
    typeBuilder.add("area", Double.class);
    SimpleFeatureType newSchema = typeBuilder.buildFeatureType();
    
    List<SimpleFeature> filteredFeatures = new ArrayList<>();
    SimpleFeatureIterator iterator = rawFeatures.features();
    
    try {
        while (iterator.hasNext()) {
            SimpleFeature feature = iterator.next();
            Geometry geom = (Geometry) feature.getDefaultGeometry();
            
            // Transform geometry to metric CRS
            Geometry projectedGeom = JTS.transform(geom, transform);
            double area = projectedGeom.getArea();
            
            // Filter by minimum area
            if (area >= minArea) {
                SimpleFeatureBuilder builder = new SimpleFeatureBuilder(newSchema);
                builder.add(projectedGeom);
                builder.add(area);
                filteredFeatures.add(builder.buildFeature(null));
            }
        }
    } finally {
        iterator.close();
    }
    
    // Export to GeoJSON
    exportToGeoJSON(filteredFeatures, newSchema);
}

private void exportToGeoJSON(List<SimpleFeature> features, SimpleFeatureType schema) throws IOException {
    FeatureJSON fjson = new FeatureJSON();
    try (StringWriter writer = new StringWriter()) {
        FeatureCollection<SimpleFeatureType, SimpleFeature> collection = new ListFeatureCollection(schema, features);
        fjson.writeFeatureCollection(collection, writer);
        System.out.println(writer.toString());
    }
}