@@ -107,13 +107,189 @@ public static Polygon toJTSPolygon(S2CellId cellId) {
107107 return new GeometryFactory ().createPolygon (coords );
108108 }
109109
110+ /**
111+ * Convert a JTS planar geometry into an S2Region whose lat/lng projection is guaranteed to
112+ * contain the input geometry.
113+ *
114+ * <p>Why a buffer is needed: Sedona geometries are planar — an edge between two vertices is a
115+ * straight line in (lng, lat) space — but S2 connects the same vertices with a great-circle arc
116+ * on the sphere. The two interpretations agree at the vertices but diverge along the edges (e.g.
117+ * the great-circle arc between two points at the same northern latitude bulges northward, leaving
118+ * the parallel that would form the planar chord). If we hand the JTS vertices to S2 directly, the
119+ * spherical polygon's interior is *smaller* than the planar polygon's interior along
120+ * non-meridional edges, so the S2 covering misses thin slivers of the original planar polygon
121+ * (see GH-2857).
122+ *
123+ * <p>To compensate, we JTS-buffer the input by an upper bound on the worst-case great-circle
124+ * deviation before converting to S2. A side effect for {@link LineString} inputs is that the
125+ * buffer turns the line into a polygon corridor; downstream callers therefore see cells in a thin
126+ * strip around the line rather than only cells the line geometrically traverses.
127+ */
110128 public static S2Region toS2Region (Geometry geom ) throws IllegalArgumentException {
111- if (geom instanceof Polygon ) {
112- return S2Utils .toS2Polygon ((Polygon ) geom );
113- } else if (geom instanceof LineString ) {
114- return S2Utils .toS2PolyLine ((LineString ) geom );
129+ if (!(geom instanceof Polygon ) && !(geom instanceof LineString )) {
130+ throw new IllegalArgumentException (
131+ "only object of Polygon, LinearRing, LineString type can be converted to S2Region" );
132+ }
133+ double eps = arcChordBufferDegrees (geom );
134+ Geometry buffered = (eps > 0 ) ? geom .buffer (eps ) : geom ;
135+ if (buffered instanceof Polygon ) {
136+ return S2Utils .toS2Polygon ((Polygon ) buffered );
137+ } else if (buffered instanceof LineString ) {
138+ // Only reachable when eps == 0 (e.g. a single-point degenerate line). Normal lines
139+ // become Polygon corridors after buffer and are handled above.
140+ return S2Utils .toS2PolyLine ((LineString ) buffered );
141+ } else if (buffered instanceof MultiPolygon && buffered .getNumGeometries () > 0 ) {
142+ // JTS buffer of self-touching geometries can collapse to MultiPolygon. We can only
143+ // hand a single S2Region back to callers, so cover the largest piece — the smaller
144+ // pieces are typically tiny artifacts of the buffer operation rather than real input.
145+ Polygon largest = (Polygon ) buffered .getGeometryN (0 );
146+ for (int i = 1 ; i < buffered .getNumGeometries (); i ++) {
147+ Polygon p = (Polygon ) buffered .getGeometryN (i );
148+ if (p .getArea () > largest .getArea ()) {
149+ largest = p ;
150+ }
151+ }
152+ return S2Utils .toS2Polygon (largest );
115153 }
116154 throw new IllegalArgumentException (
117155 "only object of Polygon, LinearRing, LineString type can be converted to S2Region" );
118156 }
157+
158+ /**
159+ * Compute the JTS buffer amount (in degrees) needed so that the spherical interpretation of the
160+ * buffered geometry fully contains the original planar geometry.
161+ *
162+ * <p>The buffer must be at least as large as the largest great-circle/chord deviation among the
163+ * edges that S2 will eventually see. Polygons and lines need different bounds because JTS buffer
164+ * affects their edges differently:
165+ *
166+ * <ul>
167+ * <li><b>Polygon</b>: each existing edge is offset perpendicularly in place; corners get
168+ * rounded into many short arcs, but no edge is dramatically lengthened. The buffered
169+ * polygon's edges therefore have ~the same length as the originals, so the original
170+ * polygon's per-edge deviation is a tight upper bound on what the buffered polygon's edges
171+ * will exhibit. We use {@link #ringMaxDeviationDegrees}.
172+ * <li><b>LineString</b>: buffering produces a corridor whose long top/bottom edges span the
173+ * line's full envelope — far longer than any individual segment when consecutive segments
174+ * are collinear (JTS often simplifies them away). Per-segment deviation severely
175+ * under-bounds the corridor's actual edge deviation. We bound by virtual edges across the
176+ * envelope via {@link #envelopeDeviationDegrees}.
177+ * </ul>
178+ *
179+ * <p>The 1.5× safety multiplier absorbs numerical error and the small additional deviation the
180+ * buffered polygon's own (slightly different) edges introduce on top of the bound.
181+ */
182+ private static double arcChordBufferDegrees (Geometry geom ) {
183+ double maxDev = 0.0 ;
184+ if (geom instanceof Polygon ) {
185+ Polygon poly = (Polygon ) geom ;
186+ maxDev = Math .max (maxDev , ringMaxDeviationDegrees (poly .getExteriorRing ().getCoordinates ()));
187+ for (int i = 0 ; i < poly .getNumInteriorRing (); i ++) {
188+ maxDev =
189+ Math .max (maxDev , ringMaxDeviationDegrees (poly .getInteriorRingN (i ).getCoordinates ()));
190+ }
191+ } else if (geom instanceof LineString ) {
192+ maxDev = envelopeDeviationDegrees (geom );
193+ }
194+ return maxDev * 1.5 ;
195+ }
196+
197+ /**
198+ * Conservative deviation upper bound for a geometry, derived from its bounding envelope rather
199+ * than its actual edges.
200+ *
201+ * <p>Used for {@link LineString} inputs because, after JTS buffer, the corridor's long edges are
202+ * not the line's segments — they connect the line's extreme endpoints (or close to it). To bound
203+ * them we probe three virtual edges across the envelope:
204+ *
205+ * <ul>
206+ * <li>The two diagonals (SW–NE and NW–SE) — diagonal great-circle arcs deviate more than
207+ * east-west arcs of the same Δλ at high latitudes, and a buffered corridor's long edges can
208+ * run in either direction depending on the line's orientation.
209+ * <li>The worst-case east-west edge at whichever latitude (top or bottom of the envelope) has
210+ * the larger absolute value — east-west arcs bulge poleward, so the deviation grows with
211+ * |latitude|, and an envelope-spanning east-west arc is what a horizontal collinear line
212+ * would buffer into.
213+ * </ul>
214+ *
215+ * <p>The max across these three bounds the deviation any corridor edge could plausibly exhibit.
216+ * This deliberately over-bounds zigzag lines whose actual corridor edges are short; the
217+ * alternative (per-segment analysis) silently fails on collinear inputs.
218+ */
219+ private static double envelopeDeviationDegrees (Geometry geom ) {
220+ Envelope env = geom .getEnvelopeInternal ();
221+ if (env .isNull ()) {
222+ return 0.0 ;
223+ }
224+ Coordinate sw = new Coordinate (env .getMinX (), env .getMinY ());
225+ Coordinate se = new Coordinate (env .getMaxX (), env .getMinY ());
226+ Coordinate ne = new Coordinate (env .getMaxX (), env .getMaxY ());
227+ Coordinate nw = new Coordinate (env .getMinX (), env .getMaxY ());
228+ // For the east-west probe, pick whichever latitude band of the envelope is further from
229+ // the equator — that's where same-Δλ great-circle arcs deviate most from the chord.
230+ double signedLat =
231+ Math .abs (env .getMinY ()) > Math .abs (env .getMaxY ()) ? env .getMinY () : env .getMaxY ();
232+ Coordinate ewWest = new Coordinate (env .getMinX (), signedLat );
233+ Coordinate ewEast = new Coordinate (env .getMaxX (), signedLat );
234+ double max = edgeDeviationDegrees (sw , ne );
235+ max = Math .max (max , edgeDeviationDegrees (nw , se ));
236+ max = Math .max (max , edgeDeviationDegrees (ewWest , ewEast ));
237+ return max ;
238+ }
239+
240+ /**
241+ * Per-edge deviation bound for a ring/path: walk consecutive vertex pairs and return the largest
242+ * single-edge great-circle/chord deviation. Used for polygon rings (exterior and interior), where
243+ * buffering preserves edge lengths and per-edge analysis is tight.
244+ */
245+ private static double ringMaxDeviationDegrees (Coordinate [] coords ) {
246+ double maxDev = 0.0 ;
247+ for (int i = 0 ; i < coords .length - 1 ; i ++) {
248+ double dev = edgeDeviationDegrees (coords [i ], coords [i + 1 ]);
249+ if (dev > maxDev ) {
250+ maxDev = dev ;
251+ }
252+ }
253+ return maxDev ;
254+ }
255+
256+ /**
257+ * Primitive deviation for a single edge: the (lng, lat) distance between the planar chord
258+ * midpoint and the great-circle arc midpoint.
259+ *
260+ * <p>Why the midpoint: a great-circle arc between two points is symmetric about its midpoint, and
261+ * the (lng, lat) deviation from the chord is maximized there. So the midpoint deviation equals
262+ * the maximum deviation along the edge — there's no need to sample multiple points.
263+ *
264+ * <p>The great-circle midpoint is computed by averaging the two endpoint S2Points (unit vectors
265+ * on the sphere) and renormalizing — a standard spherical-midpoint trick. The chord midpoint is
266+ * the plain Euclidean mean of the (lng, lat) coordinates.
267+ */
268+ private static double edgeDeviationDegrees (Coordinate a , Coordinate b ) {
269+ S2Point aPt = toS2Point (a );
270+ S2Point bPt = toS2Point (b );
271+ double mx = aPt .getX () + bPt .getX ();
272+ double my = aPt .getY () + bPt .getY ();
273+ double mz = aPt .getZ () + bPt .getZ ();
274+ double norm = Math .sqrt (mx * mx + my * my + mz * mz );
275+ if (norm < 1e-15 ) {
276+ // Antipodal endpoints — the great circle through them is not unique, so there is no
277+ // well-defined midpoint to compare against. Returning 0 effectively skips this edge;
278+ // antipodal inputs aren't realistic for S2 covering anyway.
279+ return 0.0 ;
280+ }
281+ S2LatLng midSpherical = new S2LatLng (new S2Point (mx / norm , my / norm , mz / norm ));
282+ double midSphericalLat = midSpherical .latDegrees ();
283+ double midSphericalLng = midSpherical .lngDegrees ();
284+ double midChordLat = (a .y + b .y ) / 2.0 ;
285+ double midChordLng = (a .x + b .x ) / 2.0 ;
286+ double dLat = midSphericalLat - midChordLat ;
287+ double dLng = midSphericalLng - midChordLng ;
288+ // Wrap longitude difference into [-180, 180]. Without this, an edge straddling the
289+ // antimeridian (e.g. -179° to +179°) would compute dLng ≈ 358° and produce a bogus
290+ // ~360° deviation rather than the small actual deviation.
291+ if (dLng > 180.0 ) dLng -= 360.0 ;
292+ else if (dLng < -180.0 ) dLng += 360.0 ;
293+ return Math .sqrt (dLat * dLat + dLng * dLng );
294+ }
119295}
0 commit comments