Copyright 2024 - GAEL Systems
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
Checks if the geometry pass over ±180 longitude position.
The detection of antimeridian is performed according to the distance
between longitudes positions between consecutive points of the geometry
points. The distance shall be greater than 180 to avoid revert longitude
signs around greenwich meridian 0.
It is also considered crossing antimeridian when one of the polygon
longitude is exactly to +-180 degrees.
:parameter geometry: The geometry to be controlled.
:return: True if the geometry pass over antimeridian, False otherwise.
Source code in docs/rs-server/services/common/rs_server_common/footprint_facility.py
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105 | def check_cross_antimeridian(geometry: Geometry) -> bool:
"""
Checks if the geometry pass over ±180 longitude position.
The detection of antimeridian is performed according to the distance
between longitudes positions between consecutive points of the geometry
points. The distance shall be greater than 180 to avoid revert longitude
signs around greenwich meridian 0.
It is also considered crossing antimeridian when one of the polygon
longitude is exactly to +-180 degrees.
:parameter geometry: The geometry to be controlled.
:return: True if the geometry pass over antimeridian, False otherwise.
"""
# Case of Collection of geometries (i.e. Multipolygons)
if hasattr(geometry, "geoms"):
for geom in geometry.geoms:
if check_cross_antimeridian(geom):
return True
# Path of points shall exist (Polygon or Linestring)
boundary = np.array(shapely.get_coordinates(geometry))
i = 0
while i < boundary.shape[0] - 1:
if boundary[i, 0] == 180 or boundary[i, 0] == -180 or abs(boundary[i + 1, 0] - boundary[i, 0]) > 180:
return True
i += 1
return False
|
RS-Server extension: detect only raw longitude jumps, unlike check_cross_antimeridian.
check_raw_antimeridian_jump is a local extension added to make antimeridian post-processing idempotent
and avoid reworking already split geometries. Unlike check_cross_antimeridian, this ignores coordinates
already placed on +/-180.
Source code in docs/rs-server/services/common/rs_server_common/footprint_facility.py
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134 | def check_raw_antimeridian_jump(geometry: Geometry) -> bool:
"""
RS-Server extension: detect only raw longitude jumps, unlike check_cross_antimeridian.
check_raw_antimeridian_jump is a local extension added to make antimeridian post-processing idempotent
and avoid reworking already split geometries. Unlike check_cross_antimeridian, this ignores coordinates
already placed on +/-180.
"""
# Multi geometries must be checked component by component to avoid comparing unrelated coordinate sequences.
# Unlike check_cross_antimeridian, do not fall through to shapely.get_coordinates on the whole multi geometry:
# concatenating separate components could compare the last point of one part with the first point of another
# and report a false raw longitude jump.
if hasattr(geometry, "geoms"):
for geom in geometry.geoms:
if check_raw_antimeridian_jump(geom):
return True
# No component has a raw jump; stop here to avoid comparing coordinates across independent parts.
return False
# A raw antimeridian crossing is represented by a direct longitude jump greater than 180 degrees.
boundary = np.array(shapely.get_coordinates(geometry))
i = 0
while i < boundary.shape[0] - 1:
if abs(boundary[i + 1, 0] - boundary[i, 0]) > 180:
return True
i += 1
return False
|
Elaborates linestring geometry and manages antimeridian crossing.
Source code in docs/rs-server/services/common/rs_server_common/footprint_facility.py
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503 | def rework_to_linestring_geometry(geometry: Geometry):
"""Elaborates linestring geometry and manages antimeridian crossing."""
boundaries = np.array(shapely.get_coordinates(geometry))
boundaries = np.unique(boundaries.round(decimals=1), axis=0)
if check_cross_antimeridian(geometry):
_min = min(boundaries, key=lambda point: point[0])
_max = max(boundaries, key=lambda point: point[0])
lat_at_180 = _lat_cross_antimeridian(_min, _max)
negative = [-180, lat_at_180]
positive = [180, lat_at_180]
left_antimeridian = []
right_antimeridian = []
[right_antimeridian.append(boundary) for boundary in boundaries if boundary[0] > 0]
[left_antimeridian.append(boundary) for boundary in boundaries if boundary[0] < 0]
right_antimeridian = np.concatenate((right_antimeridian, np.array([positive])), axis=0)
left_antimeridian = np.concatenate((np.array([negative]), left_antimeridian), axis=0)
reworked = shapely.multilinestrings(
[
shapely.linestrings(left_antimeridian),
shapely.linestrings(right_antimeridian),
],
)
else:
reworked = shapely.linestrings(boundaries)
return reworked
|
Rework the geometry to manage polar and antimeridian singularity.
This process implements the Polar inclusive algorithm.
The objective of this algorithm is to add the North/South Pole into
the list of coordinates of geometry polygon at the antimeridian cross.
When the geometry contains the pole the single polygon geometry including
the pole in its border point list is properly interpreted by displays
systems. When the geometry does not contain the pole, the geometry shall be
split among the antimeridian line.
:param geometry: the geometry crossing the antimeridian.
:return: the modified geometry with the closest pole included at
antimeridian crossing.
Source code in docs/rs-server/services/common/rs_server_common/footprint_facility.py
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607 | def rework_to_polygon_geometry(geometry: Geometry):
"""Rework the geometry to manage polar and antimeridian singularity.
This process implements the **Polar inclusive algorithm**.
The objective of this algorithm is to add the North/South Pole into
the list of coordinates of geometry polygon at the antimeridian cross.
When the geometry contains the pole the single polygon geometry including
the pole in its border point list is properly interpreted by displays
systems. When the geometry does not contain the pole, the geometry shall be
split among the antimeridian line.
:param geometry: the geometry crossing the antimeridian.
:return: the modified geometry with the closest pole included at
antimeridian crossing.
"""
if not check_cross_antimeridian(geometry):
return geometry
count, mixed, latitudes = _num_cross_antimeridian(geometry)
if mixed and count > 2:
logger.debug(f"WARN: Crossing antimeridian {count} times ...")
logger.debug(f" And crossing equator {_num_cross_equator(geometry)} times ...")
logger.debug("Split footprint at equator to support polar inclusion")
geometry = _split_crude_polygon_to_equator(geometry)
rwrk = []
for geom in geometry.geoms:
reworked = rework_to_polygon_geometry(geom)
if isinstance(reworked, BaseMultipartGeometry):
for geo in reworked.geoms:
rwrk.append(geo)
else:
rwrk.append(reworked)
rwrk = _merge_polygon_to_equator(rwrk)
return rwrk
# Manage case of multipolygon input
if isinstance(geometry, shapely.geometry.MultiPolygon):
rwrk = []
for geom in geometry.geoms:
reworked = rework_to_polygon_geometry(geom)
if isinstance(reworked, BaseMultipartGeometry):
raise AlreadyReworkedPolygonError("Algorithm not supported already reworked inputs.")
rwrk.append(reworked)
return shapely.geometry.MultiPolygon(rwrk)
if not isinstance(geometry, shapely.geometry.Polygon):
raise ValueError("Polygon and MultiPolygon features only are " f"supported ({type(geometry).__name__})")
boundaries = np.array(shapely.get_coordinates(geometry))
i = 0
vertical_set = False
vsign = 1
while i < boundaries.shape[0] - 1:
if abs(boundaries[i + 1, 0] - boundaries[i, 0]) > 180:
if not vertical_set:
vsign = -1 if boundaries[i, 1] < 0 else 1
vertical_set = True
hsign = -1 if boundaries[i, 0] < 0 else 1
lat = _lat_cross_antimeridian(boundaries[i], boundaries[i + 1])
boundaries = np.insert(
boundaries,
i + 1,
[
[hsign * 180, lat],
[hsign * 180, vsign * 90],
[-hsign * 180, vsign * 90],
[-hsign * 180, lat],
],
axis=0,
)
i += 5
else:
i += 1
geometry_type = type(geometry)
reworked = geometry_type(boundaries)
# When the geometry does not contain pole: Cuts the geometry among the
# antimeridian line.
if not _check_contains_pole(reworked):
reworked = _split_polygon_to_antimeridian(reworked)
else:
# Case of footprint crossing equator, antimeridian and polar zone
# Split at the equator
# Warn:
if _num_cross_equator(reworked) > 1:
reworked = _split_polygon_to_equator(reworked)
# When footprint contains overlapping, it happens at polar location.
# Polygon containing overlapping are considered invalid in shapely
# library. It includes validity check and correction methods.
# The shapely correction method extrude the overlap areas and fails
# to generate patchwork of polygons at polar area. This is probably
# due to the antimeridian crossing.
# Shapely "buffer" method fixe"s the geometry merging overlapping
# regions.
if not shapely.is_valid(reworked):
reworked = (
shapely.make_valid(reworked, method="structure", keep_collapsed=False)
if _SHAPELY_SUPPORTS_NEW_MAKE_VALID
else shapely.make_valid(reworked)
)
return shapely.buffer(reworked, 0)
|