Find polygons with in a set distance of a point

BufferedPointInPolygon.java

package com.ianturton.cookbook.filters;

import com.ianturton.cookbook.operations.BufferPoint;
import org.checkerframework.checker.units.qual.C;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.data.collection.SpatialIndexFeatureCollection;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.factory.CommonFactoryFinder;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.geotools.util.factory.GeoTools;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.opengis.filter.Filter;
import org.opengis.filter.FilterFactory2;
import org.opengis.filter.expression.Expression;
import si.uom.NonSI;
import si.uom.SI;
import systems.uom.common.USCustomary;
import tech.units.indriya.quantity.Quantities;

import javax.measure.MetricPrefix;
import javax.measure.Quantity;
import javax.measure.quantity.Length;
import java.io.File;
import java.io.IOException;

public class BufferedPointInPolygon {
	FilterFactory2 filterFactory;
	private SimpleFeatureCollection features;
	private ReferencedEnvelope env;

	public static void main(String[] args) throws IOException {
		File file = null;
		if (args.length == 0) {
			// display a data store file chooser dialog for shapefiles
			file = JFileDataStoreChooser.showOpenFile("shp", null);
			if (file == null) {
				return;
			}
		} else {
			file = new File(args[0]);
			if (!file.exists()) {
				System.err.println(file + " doesn't exist");
				return;
			}
		}

		FileDataStore store = FileDataStoreFinder.getDataStore(file);
		SimpleFeatureSource featureSource = store.getFeatureSource();
		SpatialIndexFeatureCollection idx = new SpatialIndexFeatureCollection(featureSource.getFeatures());
		GeometryFactory fac = new GeometryFactory();
		Point p = fac.createPoint(new Coordinate(22,54));
		for(int i=10;i<60;i+=10) {
			Quantity<Length> dist = Quantities.getQuantity(i, MetricPrefix.KILO(SI.METRE));
			Polygon b = (Polygon) BufferPoint.bufferPoint(dist, DefaultGeographicCRS.WGS84, p);
			System.out.println(b);
			try (SimpleFeatureIterator itr = lookup(b, idx).features()) {
				while (itr.hasNext()) {
					System.out.println(itr.next().getAttribute("sovereignt"));
				}
			}
		}
	}
	static FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
	public static SimpleFeatureCollection lookup(Geometry g, SpatialIndexFeatureCollection countries) {
		
		Filter f = ff.intersects(ff.property("the_geom"), ff.literal(g));
		SimpleFeatureCollection ret = countries.subCollection(f);


		return ret;
	}

}

BufferPoint.java

package com.ianturton.cookbook.operations;

import java.awt.geom.Point2D;
import java.util.List;

import javax.measure.Quantity;
import javax.measure.Unit;
import javax.measure.UnitConverter;
import javax.measure.quantity.Length;

import org.geotools.data.DataUtilities;
import org.geotools.feature.SchemaException;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.geotools.referencing.GeodeticCalculator;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.opengis.feature.GeometryAttribute;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.feature.type.AttributeDescriptor;
import org.opengis.feature.type.AttributeType;
import org.opengis.feature.type.GeometryType;
import org.opengis.geometry.MismatchedDimensionException;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.crs.ProjectedCRS;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;

import com.ianturton.cookbook.utilities.GenerateRandomData;

import si.uom.SI;
import systems.uom.common.USCustomary;
import tech.units.indriya.quantity.Quantities;

/**
 * Take a Point in WGS84, convert to a local projection and buffer at 1km then
 * reproject to WGS84 and print.
 * 
 * 
 * @author ian
 *
 */
public class BufferPoint {
  public static void main(String[] args) {
    SimpleFeatureType schema = null;
    SimpleFeatureType outSchema = null;
    try {
      schema = DataUtilities.createType("", "Location", "locations:Point:srid=4326," + "id:Integer" // a
                                                                                                    // number
                                                                                                    // attribute
      );
      outSchema = DataUtilities.createType("", "Location", "locations:Polygon:srid=4326," + "id:Integer" // a
                                                                                                         // number
                                                                                                         // attribute
      );
    } catch (SchemaException e) {
      // TODO Auto-generated catch block
      e.printStackTrace();
      return;
    }
    SimpleFeature feature = GenerateRandomData.createSimplePointFeature(schema);
    BufferPoint buf = new BufferPoint();
    // Quantity<Length> dist = Quantities.getQuantity(10.0,
    // MetricPrefix.KILO(SI.METRE));
    Quantity<Length> dist = Quantities.getQuantity(10.0, USCustomary.MILE);

    GeometryFactory gf = new GeometryFactory();
    /*
     * Point p = gf.createPoint(new Coordinate(-73.985460,40.748754));
     * System.out.println(p.buffer(10*0.009));
     * System.out.println(buf.bufferPoint(dist, DefaultGeographicCRS.WGS84, p));
     * System.out.println(buf.simpleBuffer(dist,p)); SimpleFeature out =
     * buf.bufferFeature(feature, dist); System.exit(1);
     */
    Point point = (Point) feature.getDefaultGeometry();
    GeodeticCalculator calc = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
    
    for (int x = -170; x < 180; x += 10) {
      for (int i = -80; i < 85; i += 10) {
        point = gf.createPoint(new Coordinate(x, i));

        calc.setStartingGeographicPoint(point.getX(), point.getY());
        calc.setDirection(0.0, dist.getValue().doubleValue());
        Point2D p2 = calc.getDestinationGeographicPoint();
        calc.setDirection(90.0, dist.getValue().doubleValue());
        Point2D p3 = calc.getDestinationGeographicPoint();

        double dy = p2.getY() - point.getY();
        double dx = p3.getX() - point.getX();
        double distance = (dy + dx) / 2.0;

        // System.out.println(dx + " " + dy + " " + distance);
        Polygon p1 = (Polygon) point.buffer(distance);

        System.out.println(buf.bufferPoint(dist, DefaultGeographicCRS.WGS84, point));
        System.out.println(p1);
      }
    }
  }

  public static SimpleFeature bufferFeature(SimpleFeature feature, Quantity<Length> distance) {
    // extract the geometry
    GeometryAttribute gProp = feature.getDefaultGeometryProperty();
    CoordinateReferenceSystem origCRS = gProp.getDescriptor().getCoordinateReferenceSystem();

    Geometry geom = (Geometry) feature.getDefaultGeometry();
    Geometry retGeom = bufferPoint(distance, origCRS, geom);
    // return a new feature containing the geom
    SimpleFeatureType schema = feature.getFeatureType();
    SimpleFeatureTypeBuilder ftBuilder = new SimpleFeatureTypeBuilder();
    ftBuilder.setCRS(origCRS);

    for (AttributeDescriptor attrib : schema.getAttributeDescriptors()) {
      AttributeType type = attrib.getType();

      if (type instanceof GeometryType) {
        String oldGeomAttrib = attrib.getLocalName();
        ftBuilder.add(oldGeomAttrib, Polygon.class);
      } else {
        ftBuilder.add(attrib);
      }
    }
    ftBuilder.setName(schema.getName());

    SimpleFeatureType nSchema = ftBuilder.buildFeatureType();
    SimpleFeatureBuilder builder = new SimpleFeatureBuilder(nSchema);
    List<Object> atts = feature.getAttributes();
    for (int i = 0; i < atts.size(); i++) {
      if (atts.get(i) instanceof Geometry) {
        atts.set(i, retGeom);
      }
    }
    return builder.buildFeature(null, atts.toArray());
  }

  /**
   * @param distance
   * @param origCRS
   * @param geom
   * @return
   */
  @SuppressWarnings("unchecked")
  public static Geometry bufferPoint(Quantity<Length> distance, CoordinateReferenceSystem origCRS, Geometry geom) {
    Geometry pGeom = geom;
    MathTransform toTransform, fromTransform = null;
    // reproject the geometry to a local projection
    Unit<Length> unit = distance.getUnit();
    if (!(origCRS instanceof ProjectedCRS)) {

      double x = geom.getCoordinate().x;
      double y = geom.getCoordinate().y;

      String code = "AUTO:42001," + x + "," + y;
      // System.out.println(code);
      CoordinateReferenceSystem auto;
      try {
        auto = CRS.decode(code);
        toTransform = CRS.findMathTransform(origCRS, auto);
        fromTransform = CRS.findMathTransform(auto, origCRS);
        pGeom = JTS.transform(geom, toTransform);
        unit = SI.METRE;
      } catch (MismatchedDimensionException | TransformException | FactoryException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }

    } else {
      unit = (Unit<Length>) origCRS.getCoordinateSystem().getAxis(0).getUnit();

    }
    UnitConverter converter = distance.getUnit().getConverterTo(unit);
    // buffer
    Geometry out = pGeom.buffer(converter.convert(distance.getValue()).doubleValue());
    Geometry retGeom = out;
    // reproject the geometry to the original projection
    if (!(origCRS instanceof ProjectedCRS)) {
      try {
        retGeom = JTS.transform(out, fromTransform);

      } catch (MismatchedDimensionException | TransformException e) {
        // TODO Auto-generated catch block
        e.printStackTrace();
      }
    }
    return retGeom;
  }
  
  public static Geometry simpleBuffer(Quantity<Length> distance,Point point) {
    GeodeticCalculator calc = new  GeodeticCalculator(DefaultGeographicCRS.WGS84);
    calc.setStartingGeographicPoint(point.getX(), point.getY());
    UnitConverter converter = distance.getUnit().getConverterTo(SI.METRE);
    double d = converter.convert(distance.getValue()).doubleValue();
    calc.setDirection(0.0, d);
    Point2D p2 = calc.getDestinationGeographicPoint();
    calc.setDirection(90.0, d);
    Point2D p3 = calc.getDestinationGeographicPoint();

    double dy = p2.getY() - point.getY();
    double dx = p3.getX() - point.getX();
    double dist = (dy + dx) / 2.0;
    return (Polygon) point.buffer(dist);
  }

  /**
   * create a buffer around the geometry, assumes the geometry is in the same
   * units as the distance variable.
   * 
   * @param geom
   *          a projected geometry.
   * @param dist
   *          a distance for the buffer in the same units as the projection.
   * @return
   */
  private Geometry buffer(Geometry geom, double dist) {

    return geom.buffer(dist);

  }
}