Examples with PreparedGeometry org.locationtech.jts.geom.prep.PreparedGeometry used on opensource projects

Example 1 with PreparedGeometry

use of org.locationtech.jts.geom.prep.PreparedGeometry in project risky by amsa-code.

the class Shapefile method load.

private Shapefile load() {
    if (state.compareAndSet(NOT_LOADED, LOADED)) {
        try {
            final List<PreparedGeometry> geometries = new ArrayList<>();
            for (String typeName : datastore.getTypeNames()) {
                SimpleFeatureSource source = datastore.getFeatureSource(typeName);
                crs = source.getBounds().getCoordinateReferenceSystem();
                final SimpleFeatureCollection features = source.getFeatures();
                SimpleFeatureIterator it = features.features();
                while (it.hasNext()) {
                    SimpleFeature feature = it.next();
                    Geometry g = (Geometry) feature.getDefaultGeometry();
                    if (bufferDistance > 0)
                        g = g.buffer(bufferDistance);
                    geometries.add(PreparedGeometryFactory.prepare(g));
                    this.features.add(feature);
                }
                it.close();
            }
            this.geometries = geometries;
        } catch (IOException e) {
            throw new RuntimeException(e);
        }
    } else if (state.get() == CLOSED)
        throw new RuntimeException("Shapefile is closed and can't be accessed");
    return this;
}

Example 2 with PreparedGeometry

use of org.locationtech.jts.geom.prep.PreparedGeometry in project risky by amsa-code.

the class ShipTypeBreakdownMain method main.

public static void main(String[] args) throws FileNotFoundException, IOException {
    // load a shapefile
    final GeometryFactory gf = new GeometryFactory();
    File file = new File("/home/dxm/temp/srr.shp");
    Map<String, Serializable> map = new HashMap<>();
    map.put("url", file.toURI().toURL());
    DataStore datastore = DataStoreFinder.getDataStore(map);
    String typeName = datastore.getTypeNames()[0];
    System.out.println(typeName);
    SimpleFeatureSource source = datastore.getFeatureSource(typeName);
    final SimpleFeatureCollection features = source.getFeatures();
    final List<PreparedGeometry> geometries = new ArrayList<>();
    SimpleFeatureIterator it = features.features();
    while (it.hasNext()) {
        SimpleFeature feature = it.next();
        Geometry g = (Geometry) feature.getDefaultGeometry();
        geometries.add(PreparedGeometryFactory.prepare(g));
    }
    // System.exit(0);
    String filename = "/media/analysis/nmea/2014/NMEA_ITU_20140815.gz";
    final Set<Long> mmsi = new HashSet<Long>();
    final Set<Long> mmsiA = new HashSet<Long>();
    final Set<Long> mmsiB = new HashSet<Long>();
    Streams.extract(Streams.nmeaFromGzip(filename)).flatMap(aisPositionsOnly).lift(Logging.<TimestampedAndLine<AisPosition>>logger().showCount().every(100000).log()).doOnNext(m -> {
        AisPosition p = m.getMessage().get().message();
        if (p.getLatitude() != null && p.getLongitude() != null && contains(gf, geometries, p.getLatitude(), p.getLongitude())) {
            long mmsiNo = m.getMessage().get().message().getMmsi();
            mmsi.add(mmsiNo);
            if (m.getMessage().get().message() instanceof AisPositionA)
                mmsiA.add(mmsiNo);
            else
                mmsiB.add(mmsiNo);
        }
    }).subscribe();
    final Map<ShipTypeClass, Set<Integer>> countsByShipType = new ConcurrentHashMap<>();
    Streams.extract(Streams.nmeaFromGzip(filename)).flatMap(aisShipStaticOnly).doOnNext(m -> {
        AisShipStatic s = m.getMessage().get().message();
        if (mmsi.contains(Long.valueOf(s.getMmsi()))) {
            boolean isClassA = s instanceof AisShipStaticA;
            ShipTypeClass shipTypeClass = new ShipTypeClass(isClassA, s.getShipType());
            if (countsByShipType.get(shipTypeClass) == null)
                countsByShipType.put(shipTypeClass, new HashSet<Integer>());
            else
                countsByShipType.get(shipTypeClass).add(s.getMmsi());
        }
    }).subscribe();
    System.out.println(countsByShipType);
    Set<String> set = new TreeSet<String>();
    int sizeA = 0;
    int sizeB = 0;
    for (Entry<ShipTypeClass, Set<Integer>> s : countsByShipType.entrySet()) {
        set.add(ShipTypeDecoder.getShipType(s.getKey().shipType) + "\t" + (s.getKey().isClassA ? "A" : "B") + "\t" + s.getValue().size());
        if (s.getKey().isClassA)
            sizeA += s.getValue().size();
        else
            sizeB += s.getValue().size();
    }
    set.stream().forEach(System.out::println);
    System.out.println("Unknown\tA\t" + (mmsiA.size() - sizeA));
    System.out.println("Unknown\tB\t" + (mmsiB.size() - sizeB));
    log.info("finished");
}

Example 3 with PreparedGeometry

use of org.locationtech.jts.geom.prep.PreparedGeometry in project qupath by qupath.

the class RoiTools method computeTiledROIs.

/**
 * Create a collection of tiled ROIs corresponding to a specified parentROI if it is larger than sizeMax, with optional overlaps.
 * <p>
 * The purpose of this is to create useful tiles whenever the exact tile size may not be essential, and overlaps may be required.
 * Tiles at the parentROI boundary will be trimmed to fit inside. If the parentROI is smaller, it is returned as is.
 *
 * @param parentROI main ROI to be tiled
 * @param sizePreferred the preferred size; in general tiles should have this size
 * @param sizeMax the maximum allowed size; occasionally it is more efficient to have a tile larger than the preferred size towards a ROI boundary to avoid creating very small tiles unnecessarily
 * @param fixedSize if true, the tile size is enforced so that complete tiles have the same size
 * @param overlap optional requested overlap between tiles
 * @return
 *
 * @see #makeTiles(ROI, int, int, boolean)
 */
public static Collection<? extends ROI> computeTiledROIs(ROI parentROI, ImmutableDimension sizePreferred, ImmutableDimension sizeMax, boolean fixedSize, int overlap) {
    ROI pathArea = parentROI != null && parentROI.isArea() ? parentROI : null;
    Rectangle2D bounds = AwtTools.getBounds2D(parentROI);
    if (pathArea == null || (bounds.getWidth() <= sizeMax.width && bounds.getHeight() <= sizeMax.height)) {
        return Collections.singletonList(parentROI);
    }
    Geometry geometry = pathArea.getGeometry();
    PreparedGeometry prepared = null;
    double xMin = bounds.getMinX();
    double yMin = bounds.getMinY();
    int nx = (int) Math.ceil(bounds.getWidth() / sizePreferred.width);
    int ny = (int) Math.ceil(bounds.getHeight() / sizePreferred.height);
    double w = fixedSize ? sizePreferred.width : (int) Math.ceil(bounds.getWidth() / nx);
    double h = fixedSize ? sizePreferred.height : (int) Math.ceil(bounds.getHeight() / ny);
    // Center the tiles
    xMin = (int) (bounds.getCenterX() - (nx * w * .5));
    yMin = (int) (bounds.getCenterY() - (ny * h * .5));
    // This can be very slow if we have an extremely large number of vertices/tiles.
    // For that reason, we try to split initially by either rows or columns if needed.
    boolean byRow = false;
    boolean byColumn = false;
    Map<Integer, Geometry> rowParents = null;
    Map<Integer, Geometry> columnParents = null;
    var envelope = geometry.getEnvelopeInternal();
    if (ny > 1 && nx > 1 && geometry.getNumPoints() > 1000) {
        // If we have a lot of points, create a prepared geometry so we can check covers/intersects quickly;
        // (for a regular geometry, it would be faster to just compute an intersection and see if it's empty)
        prepared = PreparedGeometryFactory.prepare(geometry);
        var prepared2 = prepared;
        var empty = geometry.getFactory().createEmpty(2);
        byRow = nx > ny;
        byColumn = !byRow;
        double yMin2 = yMin;
        double xMin2 = xMin;
        // Compute intersection by row so that later intersections are simplified
        if (byRow) {
            rowParents = IntStream.range(0, ny).parallel().mapToObj(yi -> yi).collect(Collectors.toMap(yi -> yi, yi -> {
                double y = yMin2 + yi * h - overlap;
                var row = GeometryTools.createRectangle(envelope.getMinX(), y, envelope.getMaxX(), h + overlap * 2);
                if (!prepared2.intersects(row))
                    return empty;
                else if (prepared2.covers(row))
                    return row;
                var temp = intersect(geometry, row);
                return temp == null ? geometry : temp;
            }));
        }
        if (byColumn) {
            columnParents = IntStream.range(0, nx).parallel().mapToObj(xi -> xi).collect(Collectors.toMap(xi -> xi, xi -> {
                double x = xMin2 + xi * w - overlap;
                var col = GeometryTools.createRectangle(x, envelope.getMinY(), w + overlap * 2, envelope.getMaxX());
                if (!prepared2.intersects(col))
                    return empty;
                else if (prepared2.covers(col))
                    return col;
                var temp = intersect(geometry, col);
                return temp == null ? geometry : temp;
            }));
        }
    }
    // Geometry local is the one we're working with for the current row or column
    // (often it's the same as the full ROI)
    Geometry geometryLocal = geometry;
    // Generate all the rectangles as geometries
    Map<Geometry, Geometry> tileGeometries = new LinkedHashMap<>();
    for (int yi = 0; yi < ny; yi++) {
        double y = yMin + yi * h - overlap;
        if (rowParents != null)
            geometryLocal = rowParents.getOrDefault(y, geometry);
        for (int xi = 0; xi < nx; xi++) {
            double x = xMin + xi * w - overlap;
            if (columnParents != null)
                geometryLocal = columnParents.getOrDefault(x, geometry);
            if (geometryLocal.isEmpty())
                continue;
            // Create the tile
            var rect = GeometryTools.createRectangle(x, y, w + overlap * 2, h + overlap * 2);
            // Use a prepared geometry if we have one to check covers/intersects & save some effort
            if (prepared != null) {
                if (!prepared.intersects(rect)) {
                    continue;
                } else if (prepared.covers(rect)) {
                    tileGeometries.put(rect, rect);
                    continue;
                }
            }
            // Checking geometryLocal.intersects(rect) first is actually much slower!
            // So add everything and filter out empty tiles later.
            tileGeometries.put(rect, geometryLocal);
        }
    }
    // Compute intersections & map to ROIs
    var plane = parentROI.getImagePlane();
    var tileROIs = tileGeometries.entrySet().parallelStream().map(entry -> intersect(entry.getKey(), entry.getValue())).filter(g -> g != null).map(g -> GeometryTools.geometryToROI(g, plane)).collect(Collectors.toList());
    // If there was an exception, the tile will be null
    if (tileROIs.size() < tileGeometries.size()) {
        logger.warn("Tiles lost during tiling: {}", tileGeometries.size() - tileROIs.size());
        logger.warn("You may be able to avoid tiling errors by calling 'Simplify shape' on any complex annotations first.");
    }
    // Remove any empty/non-area tiles
    return tileROIs.stream().filter(t -> !t.isEmpty() && t.isArea()).collect(Collectors.toList());
}

Example 4 with PreparedGeometry

use of org.locationtech.jts.geom.prep.PreparedGeometry in project risky by amsa-code.

the class Shapefile method mbr.

public Rect mbr() {
    // TODO assumes that shapefile is using WGS84?
    load();
    Rect r = null;
    for (PreparedGeometry g : geometries) {
        Coordinate[] v = g.getGeometry().getEnvelope().getCoordinates();
        System.out.println(Arrays.toString(v));
        Rect rect = new Rect(v[0].y, v[0].x, v[2].y, v[2].x);
        if (r == null)
            r = rect;
        else
            r = r.add(rect);
    }
    return r;
}

Example 5 with PreparedGeometry

use of org.locationtech.jts.geom.prep.PreparedGeometry in project risky by amsa-code.

the class ShapefileUtil method findRegionCrossingPoint.

public static TimedPosition findRegionCrossingPoint(Shapefile region, Fix fix1, Fix fix2) {
    Coordinate[] coords = new Coordinate[] { new Coordinate(fix1.lon(), fix1.lat()), new Coordinate(fix2.lon(), fix2.lat()) };
    LineString line = new GeometryFactory().createLineString(coords);
    for (PreparedGeometry g : region.geometries()) {
        if (g.crosses(line)) {
            Geometry intersection = g.getGeometry().intersection(line);
            // expecting just one point
            Coordinate coord = intersection.getCoordinate();
            double longitude = coord.x;
            double latitude = coord.y;
            Position a = Position.create(fix1.lat(), fix1.lon());
            Position b = Position.create(fix2.lat(), fix2.lon());
            Position c = Position.create(latitude, longitude);
            double ac = a.getDistanceToKm(c);
            double bc = b.getDistanceToKm(c);
            if (ac == 0) {
                return new TimedPosition(fix1.lat(), fix1.lon(), fix1.time());
            } else if (bc == 0) {
                return new TimedPosition(fix2.lat(), fix2.lon(), fix2.time());
            } else {
                // predict the timestamp based on distance from a and b
                long diff = fix2.time() - fix1.time();
                long t = Math.round(fix1.time() + ac * diff / (ac + bc));
                return new TimedPosition(latitude, longitude, t);
            }
        }
    }
    throw new RuntimeException("crossing not found");
}