Speed up Geotools Vector Grid processing

Question

I am using Geotools Vector Grid on a FeatureSource (shapefile) that been projected to "AUTO:42001". Units for the map are then in one meter increments. When the FeatureSource is a large country (like India) and the grid size is set to 10km it can take around 10 minutes to generate a grid covering the country. If I set the grid size smaller the generation time can be even larger (I tried a 1km grid and it did not complete when I left it running overnight).

Is this length of time typical for the Vector Grid, or am I doing something wrong?

Is there any way to speed up the generation of the grid?

My source data is the "ne_50m_admin_0_countries.shp" file from https://www.naturalearthdata.com/downloads/50m-cultural-vectors/50m-admin-0-countries-2/

I read the file in and select the country I want using a Query, convert to the local project and the run the GridBuilder.

File file = new File("./ne_50m_admin_0_countries.shp");
FileDataStore dataStore = FileDataStoreFinder.getDataStore(file);
SimpleFeatureSource worldMapSource = dataStore.getFeatureSource();

// Select India
Filter filter = CQL.toFilter("NAME = 'India'");
Query queryIndia = new Query(worldMapSource.getSchema().getTypeName(), filter);

FeatureCollection<SimpleFeatureType, SimpleFeature> collectionIndia = worldMapSource.getFeatures(queryIndia);
SimpleFeatureSource indiaSource = DataUtilities.source(collectionIndia);

// Project india map to local projection
SimpleFeatureSource india_LOCAL = GeoProjectionHelper.convertToLocalProjection(indiaSource);

// Set the grid size (WGS84 1.0 = 1 degree ; AUTO Local 1.0 = 1 meter) and create a bounding envelope
// that is neatly aligned with the grid size
double squareWidth = 10000.0; // 10km  (10km ~ 10min to generate)
SimpleFeatureSource indiaGrid_LOCAL = GridBuilder.getIntersectingGrid(india_LOCAL, squareWidth);

---

public class GridBuilder {

   public static SimpleFeatureSource getIntersectingGrid(SimpleFeatureSource featureSource, double squareWidth) throws IOException {
    ReferencedEnvelope gridBounds = Envelopes.expandToInclude( featureSource.getBounds(), squareWidth); 

    // Create a feature type
    SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
    tb.setName("grid");
    tb.add(
        GridFeatureBuilder.DEFAULT_GEOMETRY_ATTRIBUTE_NAME,
        Polygon.class,
        gridBounds.getCoordinateReferenceSystem());
    tb.add("id", Integer.class);
    SimpleFeatureType TYPE = tb.buildFeatureType();

    // Build the grid the custom feature builder class
    GridFeatureBuilder builder = new IntersectionBuilder(TYPE,  featureSource); 

    return Grids.createSquareGrid(gridBounds, squareWidth, -1, builder);
  }
}     

---

public class IntersectionBuilder extends GridFeatureBuilder {
  final FilterFactory2 ff2 = CommonFactoryFinder.getFilterFactory2();
  final GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();

  final SimpleFeatureSource source;
  int id = 0;

  public IntersectionBuilder(SimpleFeatureType type, SimpleFeatureSource source) {
    super(type);
    this.source = source;
  }

  public void setAttributes(GridElement el, Map<String, Object> attributes) {
    attributes.put("id", ++id);
  }

  @Override
  public boolean getCreateFeature(GridElement el) {
    // filter based on grid element geometry
    Filter filter = ff2.intersects(ff2.property("the_geom"), ff2.literal(el.toGeometry()));
    boolean result = false;

    try {
      result = !source.getFeatures(filter).isEmpty();
    } catch (IOException ex) {
      throw new IllegalStateException(ex);
    }

    return result;
  }
}

---

public class GeoProjectionHelper {

    public static SimpleFeatureSource convertToLocalProjection(SimpleFeatureSource inputFeatureSource) throws IOException {

        // extract the original CRS
        CoordinateReferenceSystem origCRS = inputFeatureSource.getFeatures().getSchema().getCoordinateReferenceSystem();

        // find centre coordinate of input
        Coordinate inputCentre = inputFeatureSource.getBounds().centre();

        // find local projection
        CoordinateReferenceSystem auto = origCRS;
        SimpleFeatureCollection outputFeatureCollection = inputFeatureSource.getFeatures();
        if (!(origCRS instanceof ProjectedCRS)) {

            String code = "AUTO:42001," + inputCentre.x + "," + inputCentre.y;

            try {
                auto = CRS.decode(code);
                System.out.println(auto.toString());

            } catch (FactoryException e) {
                e.printStackTrace();
            }

            outputFeatureCollection = new ReprojectingFeatureCollection(inputFeatureSource.getFeatures(), auto);
        }

        // wrap back to SimpleFeatureSource
        return DataUtilities.source( outputFeatureCollection );
    }

Answer 1

There are some key classes missing from your question, so this may not reflect exactly what your issue is but in my tests the majority of the execution time (5.5 sec) is writing out the shapefile of 108237 polygons.

Setup time 0.687
producing a 327.0 by 331.0 grid (108237.0)
grid time 0.662
writing time 5.562

Here's how I implemented the getIntersectingGrid class:

public static SimpleFeatureSource getIntersectingGrid(ReferencedEnvelope gridBounds, double squareWidth)
      throws IOException, NoSuchAuthorityCodeException, FactoryException, TransformException {

    double x = gridBounds.getMedian(0);
    double y = gridBounds.getMedian(1);
    // reproject to local crs, respecting axis order
    double lon;
    double lat;
    if (CRS.getAxisOrder(gridBounds.getCoordinateReferenceSystem()) == CRS.AxisOrder.EAST_NORTH) {
      lat = x;
      lon = y;
    } else {
      lat = y;
      lon = x;
    }
    String code = "AUTO:42001," + lat + "," + lon;

    CoordinateReferenceSystem auto = CRS.decode(code);
    ReferencedEnvelope bounds = gridBounds.transform(auto, true);
    bounds.expandBy(squareWidth);
    double w = Math.ceil(bounds.getWidth() / squareWidth);
    double h = Math.ceil(bounds.getHeight() / squareWidth);
    System.out.println("producing a " + w + " by " + h + " grid (" + (w * h) + ")");

    SimpleFeatureTypeBuilder tb = new SimpleFeatureTypeBuilder();
    tb.setName("grid");
    tb.add(GridFeatureBuilder.DEFAULT_GEOMETRY_ATTRIBUTE_NAME, Polygon.class, auto);
    tb.add("id", Integer.class);
    SimpleFeatureType TYPE = tb.buildFeatureType();
    GridFeatureBuilder builder = new DefaultGridFeatureBuilder(TYPE.getTypeName(), auto);

    return Grids.createSquareGrid(bounds, squareWidth, -1, builder);
  }
}

Also, if you plan to move to other countries then using AUTO:42004 may be better otherwise you can run into issues at the poles with Transverse Mercator issues.

EDIT

I think the problem is that your getCreateFeature method is doing too much work by parsing a filter and then comparing it against the features in the data source. You can go with a much simpler implementation as the target geometry doesn't change.

public static class IntersectionBuilder extends GridFeatureBuilder {
final FilterFactory2 ff2 = CommonFactoryFinder.getFilterFactory2();
final GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();

final SimpleFeature source;
Geometry geom;
int id = 0;

public IntersectionBuilder(SimpleFeatureType type, SimpleFeature source, MathTransform transform) {
  super(type);
  this.source = source;

  try {

    this.geom = JTS.transform((Geometry) source.getDefaultGeometry(), transform);
  } catch (MismatchedDimensionException | TransformException e) {
    e.printStackTrace();
  }
}



@Override
public boolean getCreateFeature(GridElement el) {
  // filter based on grid element geometry

  boolean result = el.toGeometry().intersects(geom);
  return result;
}

@Override
public void setAttributes(GridElement el, Map<String, Object> attributes) {
  attributes.put("id", ++id);

}
}

This gives me a run time of 3.98 seconds in total for India.