Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions spatial-lib/project.clj
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
(defproject nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT"
:description "A spatial library for the CMR."
:url "https://github.com/nasa/Common-Metadata-Repository/tree/master/spatial-lib"
:java-source-paths ["src/java"]
:dependencies [[nasa-cmr/cmr-common-lib "0.1.1-SNAPSHOT"]
[net.jafama/jafama "2.3.1"]
[net.mikera/core.matrix "0.54.0"]
Expand Down
2 changes: 1 addition & 1 deletion spatial-lib/src/cmr/spatial/arc.clj
Original file line number Diff line number Diff line change
Expand Up @@ -467,5 +467,5 @@


(extend-protocol d/DerivedCalculator
cmr.spatial.arc.Arc
cmr.spatial.internal.arc.Arc
(calculate-derived ^Arc [^Arc a] a))
14 changes: 7 additions & 7 deletions spatial-lib/src/cmr/spatial/arc_line_segment_intersections.clj
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,7 @@
[cmr.common.util :as u]
[cmr.spatial.line-segment :as s]
[primitive-math])
(:import cmr.spatial.arc.Arc
cmr.spatial.line_segment.LineSegment
(:import cmr.spatial.line_segment.LineSegment
cmr.spatial.point.Point))
(primitive-math/use-primitive-operators)

Expand Down Expand Up @@ -69,11 +68,11 @@

(defn line-segment-arc-intersections
"Returns a list of the points where the line segment intersects the arc."
[^LineSegment ls ^Arc arc]
[^LineSegment ls arc]

(let [ls-mbr (.mbr ls)
arc-mbr1 (.mbr1 arc)
arc-mbr2 (.mbr2 arc)
arc-mbr1 (:mbr1 arc)
arc-mbr2 (:mbr2 arc)
mbr1-intersects (m/intersects-br? :geodetic ls-mbr arc-mbr1)
mbr2-intersects (and arc-mbr2 (m/intersects-br? :geodetic ls-mbr arc-mbr2))
arc-points (a/arc->points arc)
Expand Down Expand Up @@ -132,7 +131,7 @@
"Determines if line 1 and 2 intersect. A line can be an arc or a line segment."
[line1 line2]

(if (= (type line2) Arc)
(if (= (type line2) cmr.spatial.arc.Arc)
(intersections-with-arc line1 line2)
(intersections-with-line-segment line1 line2)))

Expand All @@ -149,7 +148,8 @@
[ls point]
(s/point-on-segment? ls point))

Arc
;; Clojure Arc record type
cmr.spatial.arc.Arc
(intersections-with-arc
[arc1 arc2]
(a/intersections arc1 arc2))
Expand Down
39 changes: 17 additions & 22 deletions spatial-lib/src/cmr/spatial/cartesian_ring.clj
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,8 @@
[cmr.spatial.line-segment :as s]
[cmr.spatial.derived :as d]
[cmr.common.dev.record-pretty-printer :as record-pretty-printer])
(:import cmr.spatial.line_segment.LineSegment))
(:import
cmr.spatial.line_segment.LineSegment))
(primitive-math/use-primitive-operators)

(def external-point
Expand All @@ -30,7 +31,10 @@
line-segments

;; the minimum bounding rectangle
mbr])
mbr

;; Cached Java ring for intersection operations
java-ring])

(record-pretty-printer/enable-record-pretty-printing CartesianRing)

Expand All @@ -49,29 +53,19 @@
lines)))

(defn covers-point?
"Determines if a ring covers the given point. The algorithm works by counting the number of times
an arc between the point and a known external point crosses the ring. An even count means the point
is external. An odd count means the point is inside the ring."
[^CartesianRing ring point]

;; Only do real intersection if the mbr covers the point.
(when (mbr/cartesian-covers-point? (.mbr ring) point)
(or
;; The point is actually one of the rings points
(contains? (.point_set ring) point)
;; otherwise we'll do the real intersection algorithm
(let [;; Create the test segment
crossing-line (s/line-segment point external-point)
intersections (lines-and-line-intersections (.line_segments ring) crossing-line)]
(or (odd? (count intersections))
;; if the point itself is one of the intersections then the ring covers it
(intersections point))))))
"Determines if a ring covers the given point."
[^cmr.spatial.cartesian_ring.CartesianRing ring ^cmr.spatial.point.Point point]
;; Use cached Java ring if available, otherwise create on-the-fly
(let [java-ring (or (.java_ring ring)
(cmr.spatial.internal.ring.CartesianRing/createRing (vec (.points ring))))
java-point (cmr.spatial.shape.Point. (.lon point) (.lat point))]
(.coversPoint java-ring java-point)))
Comment thread
coderabbitai[bot] marked this conversation as resolved.

(defn ring
"Creates a new ring with the given points. If the other fields of a ring are needed. The
calculate-derived function should be used to populate it."
[points]
(->CartesianRing (mapv p/with-cartesian-equality points) nil nil nil))
(->CartesianRing (mapv p/with-cartesian-equality points) nil nil nil nil))

(defn ring->line-segments
"Determines the line-segments from the points in the ring."
Expand Down Expand Up @@ -111,5 +105,6 @@
(if (.line_segments ring)
ring
(let [^CartesianRing ring (assoc ring :point-set (set (.points ring)))
^CartesianRing ring (assoc ring :line-segments (ring->line-segments ring))]
(assoc ring :mbr (ring->mbr ring))))))
^CartesianRing ring (assoc ring :line-segments (ring->line-segments ring))
^CartesianRing ring (assoc ring :mbr (ring->mbr ring))]
(assoc ring :java-ring (cmr.spatial.internal.ring.CartesianRing/createRing (vec (.points ring))))))))
51 changes: 20 additions & 31 deletions spatial-lib/src/cmr/spatial/geodetic_ring.clj
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
[cmr.spatial.mbr :as mbr]
[cmr.spatial.point :as p]
[primitive-math])
(:import cmr.spatial.arc.Arc))
(:import))

(primitive-math/use-primitive-operators)

Expand Down Expand Up @@ -44,7 +44,9 @@
;; Three points that are not within the ring. These are used to test if a point is inside or
;; outside a ring. We generate multiple external points so that we have a backup if one external
;; point is antipodal to a point we're checking is inside a ring.
external-points])
external-points
;; Cached Java ring for intersection operations
java-ring])

(record-pretty-printer/enable-record-pretty-printing GeodeticRing)

Expand Down Expand Up @@ -81,27 +83,13 @@
" point. Ring: " (pr-str ring) " point: " (pr-str point)))))))

(defn covers-point?
"Determines if a ring covers the given point. The algorithm works by counting the number of times
an arc between the point and a known external point crosses the ring. An even count means the point
is external. An odd count means the point is inside the ring."
[ring point]
;; The pre check is necessary for rings which might contain both north and south poles
{:pre [(> (count (:external-points ring)) 0)]}

(or (and (:contains-north-pole ring) (p/is-north-pole? point))
(and (:contains-south-pole ring) (p/is-south-pole? point))
;; Only do real intersection if the mbr covers the point.
(when (mbr/geodetic-covers-point? (:mbr ring) point)
(if ((:point-set ring) point)
true ; The point is actually one of the rings points
;; otherwise we'll do the real intersection algorithm
(let [external-point (choose-external-point ring point)
;; Create the test arc
crossing-arc (a/arc point external-point)
intersections (arcs-and-arc-intersections (:arcs ring) crossing-arc)]
(or (odd-long? (count intersections))
;; if the point itself is one of the intersections then the ring covers it
(intersections point)))))))
"Determines if a ring covers the given point."
[ring ^cmr.spatial.point.Point point]
;; Use cached Java ring if available, otherwise create on-the-fly
(let [java-ring (or (:java-ring ring)
(cmr.spatial.internal.ring.GeodeticRing/createRing (vec (:points ring))))
java-point (cmr.spatial.shape.Point. (.lon point) (.lat point))]
(.coversPoint java-ring java-point)))
Comment thread
coderabbitai[bot] marked this conversation as resolved.

(defn- arcs->course-rotation-direction
"Calculates the rotation direction of the arcs of a ring. Will be one of :clockwise,
Expand All @@ -121,10 +109,10 @@
(let [courses (loop [courses (transient []) arcs arcs]
(if (empty? arcs)
(persistent! courses)
(let [^Arc a (first arcs)]
(let [a (first arcs)]
(recur (-> courses
(conj! (.initial_course a))
(conj! (.ending_course a)))
(conj! (:initial-course a))
(conj! (:ending-course a)))
(rest arcs)))))
;; Add the first turn angle on again to complete the turn
courses (conj courses (first courses))]
Expand All @@ -134,7 +122,7 @@
"Creates a new ring with the given points. If the other fields of a ring are needed. The
calculate-derived function should be used to populate it."
[points]
(->GeodeticRing (mapv p/with-geodetic-equality points) nil nil nil nil nil nil nil))
(->GeodeticRing (mapv p/with-geodetic-equality points) nil nil nil nil nil nil nil nil))

(defn contains-both-poles?
"Returns true if a ring contains both the north pole and the south pole"
Expand Down Expand Up @@ -190,12 +178,12 @@
(or (.mbr ring)
(let [arcs (ring->arcs ring)
{:keys [contains-north-pole contains-south-pole]} (ring->pole-containment ring)
br (reduce (fn [br ^Arc arc]
(let [mbr1 (.mbr1 arc)
br (reduce (fn [br arc]
(let [mbr1 (:mbr1 arc)
br (if br
(mbr/union br mbr1)
mbr1)
mbr2 (.mbr2 arc)]
mbr2 (:mbr2 arc)]
(if mbr2
(mbr/union br mbr2)
br)))
Expand Down Expand Up @@ -235,4 +223,5 @@
(assoc ring :arcs (ring->arcs ring))
(ring->pole-containment ring)
(assoc ring :mbr (ring->mbr ring))
(assoc ring :external-points (ring->external-points ring))))))
(assoc ring :external-points (ring->external-points ring))
(assoc ring :java-ring (cmr.spatial.internal.ring.GeodeticRing/createRing (vec (:points ring))))))))
57 changes: 24 additions & 33 deletions spatial-lib/src/cmr/spatial/line_string.clj
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@
(:import
(cmr.spatial.arc Arc)
(cmr.spatial.line_segment LineSegment)
(cmr.spatial.mbr Mbr)))
(cmr.spatial.mbr Mbr)
(cmr.spatial.geometry LineStringIntersections)))

(primitive-math/use-primitive-operators)

Expand Down Expand Up @@ -151,44 +152,34 @@

(defn covers-point?
"Returns true if the line covers the point"
[^LineString line point]
(let [point-set (.point_set line)
segments (.segments line)]
(or (contains? point-set point)
(u/any-true? #(point-on-segment? % point) segments))))
[^LineString line ^cmr.spatial.point.Point point]
;; Delegate to Java implementation
(let [java-line (cmr.spatial.shape.LineString.
(name (.coordinate_system line))
(vec (p/points->ords (.points line))))
java-point (cmr.spatial.shape.Point. (.lon point) (.lat point))]
(LineStringIntersections/coversPoint java-line java-point)))

(defn intersects-br?
"Returns true if the line intersects the br"
[^LineString line ^Mbr br]
(when (m/intersects-br? (.coordinate_system line) (.mbr line) br)
(if (m/single-point? br)
(covers-point? line (p/point (.west br) (.north br)))

(let [coord-sys (.coordinate_system line)]
(or
;; Does the br cover any points of the line?
(u/any-true? #(m/covers-point? coord-sys br %) (.points line))
;; Does the line contain any points of the br?
(u/any-true? #(covers-point? line %) (m/corner-points br))
;; Do any of the sides intersect?
(let [segments (.segments line)
mbr-segments (s/mbr->line-segments br)]
(loop [segments segments]
(when-let [segment (first segments)]
(let [intersects? (loop [mbr-segments mbr-segments]
(when-let [mbr-segment (first mbr-segments)]
(or (seq (asi/intersections segment mbr-segment))
(recur (rest mbr-segments)))))]
(or intersects? (recur (rest segments))))))))))))

;; Delegate to Java implementation
(let [java-line (cmr.spatial.shape.LineString.
(name (.coordinate_system line))
(vec (p/points->ords (.points line))))
java-mbr (cmr.spatial.shape.Mbr. (.west br) (.north br) (.east br) (.south br))]
(LineStringIntersections/intersectsMbr java-line java-mbr)))
(defn intersects-line-string?
"Returns true if the line string instersects the other line string"
"Returns true if the line string intersects the other line string"
[line1 line2]
(u/any-true? (fn [[s1 s2]]
(seq (asi/intersections s1 s2)))
(for [segment1 (:segments line1)
segment2 (:segments line2)]
[segment1 segment2])))
;; Delegate to Java implementation
(let [java-line1 (cmr.spatial.shape.LineString.
(name (:coordinate-system line1))
(vec (p/points->ords (:points line1))))
java-line2 (cmr.spatial.shape.LineString.
(name (:coordinate-system line2))
(vec (p/points->ords (:points line2))))]
(LineStringIntersections/intersectsLineString java-line1 java-line2)))

(extend-protocol v/SpatialValidation
cmr.spatial.line_string.LineString
Expand Down
8 changes: 7 additions & 1 deletion spatial-lib/src/cmr/spatial/lr_binary_search.clj
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,13 @@
(do
(warn "Use mbr from one of the points in the polygon because lr is not found "
(pr-str polygon))
(m/point->mbr (-> polygon :rings first :points first))))))))
;; Get first point - handle both polygon and ring
(let [first-point (if (:rings polygon)
;; It's a polygon - get first point of first ring
(-> polygon :rings first :points first)
;; It's a ring - get first point directly
(-> polygon :points first))]
(m/point->mbr first-point))))))))

(comment
(require '[cmr.spatial.kml :as kml])
Expand Down
45 changes: 15 additions & 30 deletions spatial-lib/src/cmr/spatial/mbr.clj
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,9 @@
[cmr.spatial.validation :as sv]
[cmr.spatial.messages :as msg]
[cmr.common.dev.record-pretty-printer :as record-pretty-printer])
(:import cmr.spatial.point.Point))
(:import
cmr.spatial.point.Point
[cmr.spatial.geometry MbrIntersections]))

(primitive-math/use-primitive-operators)

Expand Down Expand Up @@ -299,40 +301,23 @@

(defn non-crossing-intersects-br?
"Specialized version of intersects-br? for two mbrs that don't cross the antimeridian.
Returns true if the mbr intersects the other bounding rectangle."
[coord-sys ^Mbr m1 ^Mbr m2]
Returns true if the mbr intersects the other bounding rectangle."
[coord-sys ^cmr.spatial.mbr.Mbr m1 ^cmr.spatial.mbr.Mbr m2]
(pj/assert (not (or (crosses-antimeridian? m1)
(crosses-antimeridian? m2))))
(let [w1 (.west m1)
n1 (.north m1)
e1 (.east m1)
s1 (.south m1)
w2 (.west m2)
n2 (.north m2)
e2 (.east m2)
s2 (.south m2)
m1-touches-north? (double-approx= n1 90.0 0.0000001)
m1-touches-south? (double-approx= s1 -90.0 0.0000001)
m2-touches-north? (double-approx= n2 90.0 0.0000001)
m2-touches-south? (double-approx= s2 -90.0 0.0000001)]
(or (and (range-intersects? w1 e1 w2 e2)
(range-intersects? s1 n1 s2 n2))
(and (= coord-sys :geodetic)
(or (and m1-touches-north? m2-touches-north?)
(and m1-touches-south? m2-touches-south?))))))
;; Delegate to Java implementation
(let [java-mbr1 (cmr.spatial.shape.Mbr. (.west m1) (.north m1) (.east m1) (.south m1))
java-mbr2 (cmr.spatial.shape.Mbr. (.west m2) (.north m2) (.east m2) (.south m2))]
(MbrIntersections/nonCrossingIntersects (name coord-sys) java-mbr1 java-mbr2)))

(defn intersects-br?
"Returns true if the mbr intersects the other bounding rectangle"
[coord-sys ^Mbr mbr ^Mbr other-br]
(if (and (not (crosses-antimeridian? mbr)) (not (crosses-antimeridian? other-br)))
;; optimized case for mbrs that don't cross the antimeridian
(non-crossing-intersects-br? coord-sys mbr other-br)
(let [[m1-east m1-west] (split-across-antimeridian mbr)
[m2-east m2-west] (split-across-antimeridian other-br)]
(or (non-crossing-intersects-br? coord-sys m1-east m2-east)
(and m2-west (non-crossing-intersects-br? coord-sys m1-east m2-west))
(and m1-west (non-crossing-intersects-br? coord-sys m1-west m2-east))
(and m1-west m2-west (non-crossing-intersects-br? coord-sys m1-west m2-west))))))
[coord-sys ^cmr.spatial.mbr.Mbr mbr ^cmr.spatial.mbr.Mbr other-br]
;; Delegate to Java implementation which handles all edge cases
(let [java-mbr1 (cmr.spatial.shape.Mbr. (.west mbr) (.north mbr) (.east mbr) (.south mbr))
java-mbr2 (cmr.spatial.shape.Mbr. (.west other-br) (.north other-br)
(.east other-br) (.south other-br))]
(MbrIntersections/mbrsIntersect (name coord-sys) java-mbr1 java-mbr2)))

(defn intersections
"Returns the intersection of the two minimum bounding rectangles. This could return multiple mbrs
Expand Down
Loading