2

I'm writing a web application using GeoTools (it's Java running in a web container). I need to crop a raster file based on a shape file. The core GeoTools documentation and examples I can follow suggest that cropping is limited to Envelope (rectangular) cuts.

I came across the ImageWorker class which provides "helper methods for applying JAI operations on an image." It has a promisingly named method "mask" method ... which seems like what I want ... but it's poorly documented and seems really slow.

I'm looking for examples and/or suggestions for alternative strategies.

2

I appreciate the suggestions by Moovida and John but it turns out that GeoTools does indeed support the ability to crop a raster based on an arbitrary (....) geometry. There are no good examples out there but a colleague passed me a tip suggesting the CoverageProcessor and its "ROI" argument.

private Coverage clipImageToFeatureSource(RenderedImage image,
        ReferencedEnvelope bounds,
        FeatureSource<SimpleFeatureType, SimpleFeature> featureSource)
        throws IOException, FactoryException, MismatchedDimensionException, TransformException {
    FeatureCollection<SimpleFeatureType, SimpleFeature> collection = featureSource
            .getFeatures();

    CoordinateReferenceSystem crsFeatures = featureSource.getSchema().getCoordinateReferenceSystem();
    CoordinateReferenceSystem crsMap = bounds.getCoordinateReferenceSystem();
    boolean needsReproject = !CRS.equalsIgnoreMetadata(crsFeatures, crsMap);
    MathTransform transform = CRS.findMathTransform(crsFeatures, crsMap, true);

    FeatureIterator<SimpleFeature> iterator = collection.features();
    List<Geometry> all = new ArrayList<Geometry>();
    try {
        while (iterator.hasNext()) {
            SimpleFeature feature = iterator.next();
            Geometry geometry = (Geometry) feature.getDefaultGeometry();
            if (geometry == null)
                continue;
            if (!geometry.isSimple())
                continue;
            if (needsReproject) {
                geometry = JTS.transform(geometry, transform);
                System.out.println("Reprojected a geometry.  Result is " + geometry.toString());
            }
            Geometry intersection = geometry.intersection(JTS.toGeometry(bounds));
            if (intersection.isEmpty()) {
                continue;
            }
            //String name = (String) feature.getAttribute("NAME");
            //if (name == null)
            //   name = (String) feature.getAttribute("CNTRY_NAME");            
            if(intersection instanceof MultiPolygon) {
                MultiPolygon mp = (MultiPolygon)intersection;
                for (int i = 0; i < mp.getNumGeometries(); i++) {
                    com.vividsolutions.jts.geom.Polygon g = (com.vividsolutions.jts.geom.Polygon)mp.getGeometryN(i);
                    Geometry gIntersection = IntersectUtils.intersection(g, JTS.toGeometry(bounds));
                    if (gIntersection.isEmpty()) {
                        continue;
                    }
                    all.add(g);
                }
            }
            else if (intersection instanceof Polygon)
                all.add(intersection);
            else
                continue;
        }
    } finally {
        if (iterator != null) {
            iterator.close();
        }
    }
    GridCoverageFactory gridCoverageFactory = new GridCoverageFactory();
    Coverage coverage = gridCoverageFactory.create("Raster", image, bounds);
    Coverage clippedCoverage = null;
    if (all.size() > 0) {
        CoverageProcessor processor = new CoverageProcessor();
        ParameterValueGroup params = processor.getOperation("CoverageCrop")
                .getParameters();
        params.parameter("Source").setValue(coverage);
        GeometryFactory factory = JTSFactoryFinder.getGeometryFactory(null);
        Geometry[] a = all.toArray(new Geometry[0]);
        GeometryCollection c = new GeometryCollection(a, factory);
        //params.parameter("ENVELOPE").setValue(bounds);
        params.parameter("ROI").setValue(c);
        params.parameter("ForceMosaic").setValue(true);
        clippedCoverage = processor.doOperation(params);
    }
    if (all.size() == 0){
        logger.info("Crop by shapefile requested but no simple features matched extent!");
    }
    return clippedCoverage;
}

This clips the image but may also reduce the original extent. If you want to preserve that, then you'll need to "matt" the clippedCoverage like this:

    private BufferedImage mattCroppedImage(final BufferedImage source, GridCoverage2D cropped) 
    {
    RenderedImage raster = cropped.getRenderedImage();
    int height = source.getHeight();
            int width = source.getWidth();
            BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_INT_RGB);
            Graphics2D gr = image.createGraphics();
            gr.setPaint(Color.green);
    gr.fill(new Rectangle2D.Double(0,0, image.getWidth(), image.getHeight()));
    AffineTransform at = AffineTransform.getTranslateInstance(0, 0);
    gr.drawRenderedImage(cropped.getRenderedImage(), at);
    return image;
}
  • Now that I remember, the ZonalStats operation in geotools does exactly that to get info from shapes. I think you should be able to get what you need from the testcase. In there you have also a method to get to an AWT polygon from the "geographic" one, which should then be the ROI. – moovida Oct 5 '13 at 20:53
  • Edited to show more robust solution. – msr Nov 21 '13 at 19:39
  • I used your approach, in my case however, it clipped to the bounding box of the shape, not the shape itself. I had to do a refined test, i.e., I had to check each pixel if its in the polygon (after clipping). – Mr M Apr 24 '16 at 11:43
0

I know you've stated that you are using GeoTools, but if you're willing to use other software as an alternative strategy, then I would suggest using the 'Clip Raster to Polygon' tool in the open-source GIS Whitebox GAT (http://www.uoguelph.ca/~hydrogeo/Whitebox/). It does exactly that, crops a raster to a vector polygon file. Full disclosure: I am the lead developer of Whitebox.

enter image description here

0

If the image needs to be cropped to any shape, you might want to use a rasterization and masking approach.

One example of doing this with geotools is done in the jgrasstools ScanlineRasterizer:

https://code.google.com/p/jgrasstools/source/browse/jgrassgears/src/main/java/org/jgrasstools/gears/modules/r/scanline/OmsScanLineRasterizer.java

The result is a raster mask, that can be used to extract only the needed part.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.