crs-l2-fiji-antimeridian

Rerun

import geopandas as gpd

# Load the data
gdf = gpd.read_file('fiji_transects_wgs84.geojson')

# Inspect the data
print("Columns:", gdf.columns.tolist())
print("Geometry types:", gdf.geometry.type.value_counts())
print("First few rows:")
print(gdf.head())
print("\nCRS:", gdf.crs)

import geopandas as gpd
import pandas as pd
from shapely.geometry import LineString, MultiLineString
from shapely.ops import transform
import pyproj

def fix_antimeridian_geometry(line, lon_range=(-180, 180)):
    """
    If a LineString crosses the antimeridian, split it into multiple LineStrings.
    This function assumes the geometry might have long-distance jumps if it was 
    stored incorrectly or if it crosses the 180/-180 boundary.
    However, the user says they 'render as arcs spanning the entire globe'.
    This usually happens when a line goes from e.g. 179 to -179 and is 
    interpreted as a great circle or a long segment.
    """
    # Check if the line is a single segment that jumps across the antimeridian
    # We'll look at the coordinates.
    coords = list(line.coords)
    
    # If the difference in longitude between consecutive points is huge, 
    # it's likely a jump across the antimeridian.
    # But the user says they "render as arcs spanning the entire globe".
    # This often means the coordinates are 'wrapped' incorrectly or the 
    # projection used when creating them was bad.
    
    # In GeoJSON, if we have a line from 179 to -179, it's a long segment.
    # A proper way to fix this for a local area is to shift the coordinates.
    # However, we don't know the center. 
    # Let's look at the bounding box.
    
    return line

def split_antimeridian_lines(gdf):
    """
    Splits lines that cross the antimeridian into multiple parts.
    For Fiji, they are near the antimeridian (around 178-180).
    We can shift everything so they are centered around 0 if necessary, 
    but better yet, we can detect if the longitude jump is > 180.
    """
    new_geoms = []
    new_indices = []
    
    for idx, row in gdf.iterrows():
        geom = row.geometry
        if geom.geom_type == 'LineString':
            coords = list(geom.coords)
            parts = []
            current_part = [coords[0]]
            
            for i in range(1, len(coords)):
                prev_lon = coords[i-1][0]
                curr_lon = coords[i][0]
                
                # Detect jump
                diff = curr_lon - prev_lon
                if diff > 180:
                    # Jump from positive to negative
                    # The segment is actually short, crossing 180
                    # We can 'wrap' the curr_lon
                    actual_curr_lon = curr_lon - 360
                    # Wait, if we just want to split it into two lines:
                    # The segment from prev to curr is the problem.
                    # But if it's a single segment, we split it at the boundary.
                    # However, the instruction says "the problematic ones must end up as multi-part geometries".
                    # So we split the line into two segments: [prev, 180] and [-180, curr]
                    # But the geometry itself might be a single segment in the file.
                    
                    # Let's try a simpler approach:
                    # If a line crosses the antimeridian, we split it.
                    pass
            # Re-thinking: 
            # If the user says "render as arcs spanning the entire globe", 
            # it's because the longitude jumps from ~180 to ~-180.
            # Let's check for jumps.
            
            coords = list(geom.coords)
            is_crossing = False
            for i in range(len(coords)-1):
                if abs(coords[i][0] - coords[i+1][0]) > 180:
                    is_crossing = True
                    break
            
            if is_crossing:
                # Split into multiple LineStrings
                # We'll use a simple heuristic: split at the jump.
                # For a single jump, we'll create two segments.
                # But we need to be careful about the coordinate values.
                # If it goes from 179 to -179, we split it into 179->180 and -180->-179?
                # No, actually we should split the geometry at the point where it crosses.
                # But we don't have the point. 
                # A better way: if the jump is > 180, it's a jump.
                
                # Let's use a more robust approach:
                # If we treat it as a "wrapped" geometry, we can use a projection that handles it,
                # or we can manually split the segments.
                
                # Let's try to split the LineString at the jump.
                new_parts = []
                temp_coords = []
                for i in range(len(coords)):
                    if i > 0:
                        prev_lon = coords[i-1][0]
                        curr_lon = coords[i][0]
                        if abs(prev_lon - curr_lon) > 180:
                            # Jump detected!
                            # Split the current part and start a new one.
                            if temp_coords:
                                new_parts.append(LineString(temp_coords))
                                temp_coords = []
                            
                            # For the jump segment:
                            # The segment itself is what causes the "arc".
                            # But a LineString can only have one segment between two points.
                            # If the segment is (179, -179), it's one segment.
                            # To "fix" it, we must split it into two segments? 
                            # But a LineString can't have a "gap" unless it's a MultiLineString.
                            # So we'll create two segments: (179, 180) and (-180, -179) 
                            # but we don't know the exact crossing.
                            # Actually, if we're in Fiji, they are likely near 180.
                            
                            # Wait, if the segment is (179, -179), we can't just split it 
                            # into two segments without adding a point at 180/-180.
                            
                            # Let's try a different approach:
                            # 1. Shift all longitudes to be within [0, 360]
                            # 2. Then check for jumps.
                            # 3. Or even better: use a projection that's local.
                            pass

    return gdf

# Let's try a simpler logic for the specific problem described.
# The user says they render as arcs. This is because of the jump.
# If I split the LineString at the jump and "wrap" the coordinates.

def fix_geometry_v2(geom):
    if geom.geom_type != 'LineString':
        return geom
    
    coords = list(geom.coords)
    new_coords = []
    
    # To handle the wrap-around, we'll check if the jump is > 180.
    # We'll transform the coordinates to a continuous range.
    # e.g. if it goes 179 -> -179, we treat it as 179 -> 181.
    
    continuous_coords = [coords[0]]
    for i in range(1, len(coords)):
        prev_lon, prev_lat = continuous_coords[-1]
        curr_lon, curr_lat = coords[i]
        
        diff = curr_lon - prev_lon
        # If jump is > 180, it's a wrap from 180 to -180
        if diff > 180:
            diff -= 360
        # If jump is < -180, it's a wrap from -180 to 180
        elif diff < -180:
            diff += 360
            
        continuous_coords.append((prev_lon + diff, curr_lat))
    
    # Now continuous_coords are in a continuous longitudinal space.
    # We need to split this into multiple LineStrings if it crosses the 180 boundary.
    # However, the problem is we need to return them as MultiLineString in the 
    # original -180 to 180 space.
    
    # Step 1: Detect where it crosses the 180 / -180 boundary in continuous space.
    # Actually, we can just work with the continuous coords and then 
    # wrap them back to -180, 180.
    
    # Let's find the boundaries.
    # If we have continuous coords like [(179, -17), (181, -17), (182, -17)]
    # We should split it into [(179, -17), (180, -17)] and [(-180, -17), (-178, -17)]
    
    # A more robust way:
    # 1. Convert to continuous longitudes.
    # 2. For each segment (p1, p2) in continuous longitudes:
    #    If the segment crosses a multiple of 180:
    #       Split it.
    # 3. Wrap all resulting segments back to [-180, 180].
    
    # Let's implement this.
    pass

def robust_fix(geom):
    if geom.geom_type != 'LineString':
        return geom
    
    coords = list(geom.coords)
    if not coords:
        return geom
        
    # 1. Create continuous coordinates
    cont_coords = [coords[0]]
    for i in range(1, len(coords)):
        prev_lon, prev_lat = cont_coords[-1]
        curr_lon, curr_lat = coords[i]
        
        diff = curr_lon - prev_lon
        if diff > 180:
            diff -= 360
        elif diff < -180:
            diff += 360
        
        cont_coords.append((prev_lon + diff, curr_lat))
        
    # 2. Split continuous coords into segments that don't cross the 180 boundary
    # The boundary is at ... -180, 180, 540, ...
    # Since we are in Fiji, we only care about 180 (or -180).
    
    lines = []
    current_segment = [cont_coords[0]]
    
    for i in range(1, len(cont_coords)):
        p1 = cont_coords[i-1]
        p2 = cont_coords[i]
        
        # Check if segment [p1, p2] crosses any multiple of 180? 
        # Actually, it's simpler: just check if the longitude crosses the 180 mark.
        # Because we used continuous longitudes, the 180 boundary is at 180, 540, -180, etc.
        # A segment crosses a boundary if floor(lon/180) changes? 
        # No, the boundary is specifically at 180, 540, -180, etc.
        # But in -180 to 180, the only "boundary" is 180 or -180.
        # In continuous space, it's 180 + 360*k.
        
        # Let's check if there's an integer 'k' such that 180 + 360*k is between p1[0] and p2[0].
        # Since our jumps are <= 180, there can be at most one such k.
        
        # A simpler way to check for crossing:
        # Does the segment cross (180 + 360*k) for any k?
        # Since we're looking at the interval [p1[0], p2[0]], we can check:
        # Is there an integer k such that p1[0] <= 180 + 360*k <= p2[0] (or vice versa)?
        
        cross_boundary = False
        # Possible k's:
        # k = round((p1[0] - 180) / 360)
        # Let's check k-1, k, k+1
        for k in range(-2, 3): # Sufficient for small segments
            boundary = 180 + 360 * k
            if (p1[0] < boundary <= p2[0]) or (p2[0] < boundary <= p1[0]):
                # Crossing detected!
                # We need to split the segment at the boundary.
                # We can't just add the boundary point to the LineString if it's not in the original.
                # Wait, the instruction says "the problematic ones must end up as multi-part geometries that faithfully represent the actual survey path."
                # If the line is (179, -17) to (-179, -17), it's ONE segment in the file.
                # To "faithfully represent" it, we should split it at the 180/-180 line.
                
                # Segment is (p1_lon, p1_lat) to (p2_lon, p2_lat) in continuous space.
                # Boundary is 'boundary'.
                # Interpolate the latitude at the boundary.
                ratio = (boundary - p1[0]) / (p2[0] - p1[0])
                boundary_lat = p1[1] + ratio * (p2[1] - p1[1])
                
                # Split the segment into (p1, boundary_point) and (boundary_point, p2)
                # Actually, we'll just finish the current segment with p1, and start a new one with p2.
                # But that's not quite right for a single segment.
                # If a single segment is (179, -17) to (-179, -17), 
                # it should be split into (179, -17) -> (180, -17) and (-180, -17) -> (-179, -17).
                # But we don't know the original segment's latitude at 180. 
                # We interpolate it.
                
                # Let's refine:
                # We split the continuous line into segments at every 180 + 360k.
                # For each continuous segment (p1, p2) that crosses a boundary:
                #   Create segment (p1, boundary_point)
                #   Create segment (boundary_point, p2)
                # This might add points not in the original, but it's necessary to avoid the "arc".
                
                # Actually, if it's a single LineString segment in the input, we'll just split it.
                pass

    # Let's simplify. We can use a library or a more direct method.
    # Let's use the fact that we're in Fiji. Most coords are around 175-180 or -175 to -180.
    
    # Let's try this:
    # 1. Get continuous coords.
    # 2. For each segment in continuous coords:
    #    If it crosses a boundary:
    #       Split it into two parts: (start_of_segment, boundary) and (boundary, end_of_segment).
    #    Else:
    #       Keep it as one segment.
    # 3. Collect all resulting segments.
    # 4. Wrap each segment's coords back to [-180, 180].
    # 5. Group them into MultiLineStrings.
    
    # Wait, the problem says "the problematic ones must end up as multi-part geometries".
    # This means a LineString that was problematic becomes a MultiLineString.
    
    # Let's write the robust logic.
    pass

def solve():
    import geopandas as gpd
    from shapely.geometry import LineString, MultiLineString
    import numpy as np

    gdf = gpd.read_file('fiji_transects_wgs84.geojson')
    
    new_geometries = []
    
    for idx, row in gdf.iterrows():
        geom = row.geometry
        if geom.geom_type != 'LineString':
            new_geometries.append(geom)
            continue
        
        coords = list(geom.coords)
        if len(coords) < 2:
            new_geometries.append(geom)
            continue
            
        # 1. Continuous coordinates
        cont_coords = [coords[0]]
        for i in range(1, len(coords)):
            prev_lon, prev_lat = cont_coords[-1]
            curr_lon, curr_lat = coords[i]
            diff = curr_lon - prev_lon
            if diff > 180:
                diff -= 360
            elif diff < -180:
                diff += 360
            cont_coords.append((prev_lon + diff, curr_lat))
            
        # 2. Split at boundaries
        segments = []
        current_seg = [cont_coords[0]]
        
        for i in range(1, len(cont_coords)):
            p1 = cont_coords[i-1]
            p2 = cont_coords[i]
            
            # Check for boundary crossing in this segment
            # Boundary is 180 + 360*k
            # Since it's a single segment, we check if any boundary is between p1[0] and p2[0]
            # We can find k by checking the range of possible k's.
            # But it's even simpler: a boundary is crossed if floor(p1/360) != floor(p2/360) 
            # No, the boundary is at 180 + 360k.
            # Let's check if there's an integer k such that p1_lon < 180 + 360k <= p2_lon (if p2 > p1)
            
            found_boundary = False
            # The boundaries are 180, 540, -180, -540...
            # We can just find if any (180 + 360*k) is in (min(p1,p2), max(p1,p2)]
            low, high = min(p1[0], p2[0]), max(p1[0], p2[0])
            
            # Solve for k: low < 180 + 360k <= high
            # (low - 180) / 360 < k <= (high - 180) / 360
            k_min = int(np.floor((low - 180) / 360))
            k_max = int(np.floor((high - 180) / 360))
            
            for k in range(k_min - 1, k_max + 2):
                boundary = 180 + 360 * k
                if low < boundary <= high:
                    found_boundary = True
                    # Interpolate
                    ratio = (boundary - p1[0]) / (p2[0] - p1[0])
                    boundary_lat = p1[1] + ratio * (p2[1] - p1[1])
                    boundary_pt = (boundary, boundary_lat)
                    
                    # The segment is split into (p1, boundary_pt) and (boundary_pt, p2)
                    # But wait, 'current_seg' might have multiple points.
                    # We should close the current segment with 'boundary_pt' (not exactly, but 
                    # the point where the boundary was crossed).
                    
                    # Actually, if the crossing happens *inside* the segment (p1, p2):
                    # We finish the current_seg with p1, then we need a new segment starting from 
                    # a point 'at the boundary'. 
                    # However, if p1 is the start of a segment and p2 is the end, and they cross...
                    # We should add the boundary point to the current segment.
                    
                    # Let's handle this more cleanly.
                    pass
            
            # Revised strategy for splitting:
            # Instead of segment by segment, let's find all boundary points 
            # in the continuous sequence and insert them.
            
        # Let's do the insertion approach.
        all_cont_coords = []
        for i in range(len(cont_coords) - 1):
            p1 = cont_coords[i]
            p2 = cont_coords[i+1]
            all_cont_coords.append(p1)
            
            # Check for boundary in (p1, p2]
            low, high = min(p1[0], p2[0]), max(p1[0], p2[0])
            k_min = int(np.floor((low - 180) / 360))
            k_max = int(np.floor((high - 180) / 360))
            for k in range(k_min - 1, k_max + 2):
                boundary = 180 + 360 * k
                if low < boundary <= high:
                    ratio = (boundary - p1[0]) / (p2[0] - p1[0])
                    boundary_lat = p1[1] + ratio * (p2[1] - p1[1])
                    all_cont_coords.append((boundary, boundary_lat))
                    
        all_cont_coords.append(cont_coords[-1])
        
        # Now we have a continuous sequence of points, some of which are on the boundary.
        # We split this sequence into multiple LineStrings whenever we cross a boundary.
        # A boundary is a point where the longitude is 180 + 360k.
        
        split_segments = []
        if all_cont_coords:
            current_seg_coords = [all_cont_coords[0]]
            for i in range(1, len(all_cont_coords)):
                prev_pt = all_cont_coords[i-1]
                curr_pt = all_cont_coords[i]
                
                # If prev_pt is a boundary, we start a new segment.
                # Note: A boundary point can be the end of one segment and the start of another.
                # To avoid zero-length segments, we check.
                is_boundary = False
                # Use a small epsilon for float comparison
                if abs((prev_pt[0] - 180) % 360) < 1e-9 or abs((prev_pt[0] + 180) % 360) < 1e-9:
                    # It's a boundary. 
                    # However, the boundary point is also the start of the NEXT segment.
                    # If it's a boundary point, we should start a new segment.
                    # But if it's the very first point, we don't split.
                    if i > 1: # Not the first point
                        # Check if it's actually a boundary (180 + 360k)
                        # We'll use a more robust check: is it a multiple of 180?
                        # Wait, the boundaries are ... -180, 180, 540 ...
                        # If it's 180, it's a boundary. If it's 540, it's a boundary.
                        # If it's 0, it's not.
                        # So (pt[0] - 180) % 360 == 0.
                        if abs((prev_pt[0] - 180) % 360) < 1e-7:
                            is_boundary = True
                
                if is_boundary:
                    # Close current segment and start new one
                    if len(current_seg_coords) >= 2:
                        split_segments.append(LineString(current_seg_coords))
                    current_seg_coords = [curr_pt]
                else:
                    current_seg_coords.append(curr_pt)
            
            if len(current_seg_coords) >= 2:
                split_segments.append(LineString(current_seg_coords))

        # 3. Wrap back to [-180, 180]
        final_segments = []
        for seg in split_segments:
            wrapped_coords = []
            for lon, lat in seg.coords:
                # Wrap lon to [-180, 180]
                # (lon + 180) % 360 - 180
                w_lon = ((lon + 180) % 360 + 360) % 360 - 180
                # Handle the -180 vs 180 case: if w_lon is -180, it's fine. 
                # But if it was 180, it's also fine.
                wrapped_coords.append((w_lon, lat))
            
            if len(wrapped_coords) >= 2:
                final_segments.append(LineString(wrapped_coords))
        
        if final_segments:
            if len(final_segments) == 1:
                new_geometries.append(final_segments[0])
            else:
                new_geometries.append(MultiLineString(final_segments))
        else:
            # If for some reason we have no segments, just use original (though shouldn't happen)
            new_geometries.append(geom)

    # Reconstruct GDF
    new_gdf = gdf.copy()
    new_gdf.geometry = new_geometries
    return new_gdf

# Actually, the requirement is that EVERY feature must be a MultiLineString.
def force_multilinestring(gdf):
    new_geoms = []
    for geom in gdf.geometry:
        if geom.geom_type == 'LineString':
            new_geoms.append(MultiLineString([geom]))
        else:
            new_geoms.append(geom)
    gdf.geometry = new_geoms
    return gdf

# Now add length_m
def add_length(gdf):
    # Convert to Fiji's national metric grid first.
    # The user said "convert them to Fiji's national metric grid".
    # Fiji's national grid is usually a UTM zone.
    # Fiji is in UTM zones 58S and 59S. 
    # However, "Fiji's national metric grid" might refer to a specific CRS.
    # If not specified, I'll use a suitable UTM zone for Fiji or a projected CRS.
    # Let's check if there's a standard. Common is EPSG:32758 (UTM 58S) or 32759 (UTM 59S).
    # A more general choice is a local projection.
    # Let's use EPSG:32758 (UTM zone 58S) as it covers a large part of Fiji.
    # Or better, I'll check if the data is mostly in one zone.
    
    # Let's assume EPSG:32758 for now, or use a more general one if needed.
    # Actually, many people use EPSG:3857 for web, but that's not metric for Fiji.
    # Let's try to find the best UTM zone.
    
    # To be safe, I'll project to a metric CRS.
    # Fiji is around -18, 178.
    # EPSG:32758 is UTM 58S (covers 174E-180E).
    # EPSG:32759 is UTM 59S (covers 180E-186E).
    # Since the data is around 178-179, 58S is appropriate.
    
    gdf_metric = gdf.to_crs(epsg=32758)
    gdf_metric['length_m'] = gdf_metric.geometry.length
    return gdf_metric

# Wait, the user wants the result as GeoJSON named fiji_transects_fmg.geojson.
# And the length should be in the GeoJSON.
# The final output should be in Fiji's national metric grid.
# If it's in a metric grid, the CRS will be something like EPSG:32758.

# Let's check the CRS of the input. It's EPSG:4326.
# The result should be in the metric grid.
# So I'll:
# 1. Fix the antimeridian geometries (stay in 4326 for now).
# 2. Project to Fiji's metric grid (e.g., EPSG:32758).
# 3. Calculate length.
# 4. Force MultiLineString.
# 5. Save.

# Let's refine the split logic one more time.
# If I have a LineString (179, -17) -> (-179, -17).
# Continuous: (179, -17) -> (181, -17).
# Boundary is at 180.
# Points: (179, -17), (180, -17), (181, -17).
# Segments: [(179, -17), (180, -17)] and [(180, -17), (181, -17)].
# Wrapped: [(179, -17), (180, -17)] and [(-180, -17), (-179, -17)].
# This looks correct.

# Let's execute.
import geopandas as gpd
import numpy as np
from shapely.geometry import LineString, MultiLineString

def fix_antimeridian(gdf):
    new_geoms = []
    for idx, row in gdf.iterrows():
        geom = row.geometry
        if geom.geom_type != 'LineString':
            new_geoms.append(geom)
            continue
        
        coords = list(geom.coords)
        if len(coords) < 2:
            new_geoms.append(geom)
            continue
            
        cont_coords = [coords[0]]
        for i in range(1, len(coords)):
            prev_lon, prev_lat = cont_coords[-1]
            curr_lon, curr_lat = coords[i]
            diff = curr_lon - prev_lon
            if diff > 180: diff -= 360
            elif diff < -180: diff += 360
            cont_coords.append((prev_lon + diff, curr_lat))
            
        all_cont_coords = []
        for i in range(len(cont_coords) - 1):
            p1 = cont_coords[i]
            p2 = cont_coords[i+1]
            all_cont_coords.append(p1)
            low, high = min(p1[0], p2[0]), max(p1[0], p2[0])
            k_min = int(np.floor((low - 180) / 360))
            k_max = int(np.floor((high - 180) / 360))
            for k in range(k_min - 1, k_max + 2):
                boundary = 180 + 360 * k
                if low < boundary <= high:
                    ratio = (boundary - p1[0]) / (p2[0] - p1[0])
                    boundary_lat = p1[1] + ratio * (p2[1] - p1[1])
                    all_cont_coords.append((boundary, boundary_lat))
        all_cont_coords.append(cont_coords[-1])
        
        split_segments = []
        if all_cont_coords:
            current_seg_coords = [all_cont_coords[0]]
            for i in range(1, len(all_cont_coords)):
                prev_pt = all_cont_coords[i-1]
                curr_pt = all_cont_coords[i]
                # Boundary check: is prev_pt a boundary?
                is_boundary = False
                if i > 1:
                    # If (prev_pt[0] - 180) % 360 is close to 0
                    if abs((prev_pt[0] - 180) % 360) < 1e-7 or abs((prev_pt[0] + 180) % 360) < 1e-7:
                        # But we need to be sure it's an actual 180, 540 etc.
                        # Actually, (prev_pt[0] - 180) % 360 is enough.
                        if abs((prev_pt[0] - 180) % 360) < 1e-7:
                            is_boundary = True
                
                if is_boundary:
                    if len(current_seg_coords) >= 2:
                        split_segments.append(LineString(current_seg_coords))
                    current_seg_coords = [curr_pt]
                else:
                    current_seg_coords.append(curr_pt)
            if len(current_seg_coords) >= 2:
                split_segments.append(LineString(current_seg_coords))
        
        final_segments = []
        for seg in split_segments:
            wrapped_coords = []
            for lon, lat in seg.coords:
                w_lon = ((lon + 180) % 360 + 360) % 360 - 180
                # Handle the 180 boundary case carefully: 
                # if w_lon is -180 but the original was 180, it doesn't matter much for geometry.
                wrapped_coords.append((w_lon, lat))
            if len(wrapped_coords) >= 2:
                final_segments.append(LineString(wrapped_coords))
        
        if final_segments:
            if len(final_segments) == 1:
                new_geoms.append(final_segments[0])
            else:
                new_geoms.append(MultiLineString(final_segments))
        else:
            new_geoms.append(geom)
            
    gdf.geometry = new_geoms
    return gdf

# After running this, I'll wrap it in the final script.

import geopandas as gpd
import numpy as np
from shapely.geometry import LineString, MultiLineString

def fix_antimeridian(gdf):
    new_geoms = []
    for idx, row in gdf.iterrows():
        geom = row.geometry
        if geom.geom_type != 'LineString':
            new_geoms.append(geom)
            continue
        
        coords = list(geom.coords)
        if len(coords) < 2:
            new_geoms.append(geom)
            continue
            
        cont_coords = [coords[0]]
        for i in range(1, len(coords)):
            prev_lon, prev_lat = cont_coords[-1]
            curr_lon, curr_lat = coords[i]
            diff = curr_lon - prev_lon
            if diff > 180: diff -= 360
            elif diff < -180: diff += 360
            cont_coords.append((prev_lon + diff, curr_lat))
            
        all_cont_coords = []
        for i in range(len(cont_coords) - 1):
            p1 = cont_coords[i]
            p2 = cont_coords[i+1]
            all_cont_coords.append(p1)
            low, high = min(p1[0], p2[0]), max(p1[0], p2[0])
            k_min = int(np.floor((low - 180) / 360))
            k_max = int(np.floor((high - 180) / 360))
            for k in range(k_min - 1, k_max + 2):
                boundary = 180 + 360 * k
                if low < boundary <= high:
                    ratio = (boundary - p1[0]) / (p2[0] - p1[0])
                    boundary_lat = p1[1] + ratio * (p2[1] - p1[1])
                    all_cont_coords.append((boundary, boundary_lat))
        all_cont_coords.append(cont_coords[-1])
        
        split_segments = []
        if all_cont_coords:
            current_seg_coords = [all_cont_coords[0]]
            for i in range(1, len(all_cont_coords)):
                prev_pt = all_cont_coords[i-1]
                curr_pt = all_cont_coords[i]
                is_boundary = False
                if i > 1:
                    if abs((prev_pt[0] - 180) % 360) < 1e-7:
                        is_boundary = True
                
                if is_boundary:
                    if len(current_seg_coords) >= 2:
                        split_segments.append(LineString(current_seg_coords))
                    current_seg_coords = [curr_pt]
                else:
                    current_seg_coords.append(curr_pt)
            if len(current_seg_coords) >= 2:
                split_segments.append(LineString(current_seg_coords))
        
        final_segments = []
        for seg in split_segments:
            wrapped_coords = []
            for lon, lat in seg.coords:
                w_lon = ((lon + 180) % 360 + 360) % 360 - 180
                wrapped_coords.append((w_lon, lat))
            if len(wrapped_coords) >= 2:
                final_segments.append(LineString(wrapped_coords))
        
        if final_segments:
            if len(final_segments) == 1:
                # The requirement says EVERY feature must be MultiLineString.
                # So even if there is one segment, it must be wrapped in MultiLineString.
                new_geoms.append(MultiLineString([final_segments[0]]))
            else:
                new_geoms.append(MultiLineString(final_segments))
        else:
            new_geoms.append(geom)
            
    gdf.geometry = new_geoms
    return gdf

def main():
    # Load
    gdf = gpd.read_file('fiji_transects_wgs84.geojson')
    
    # 1. Fix antimeridian
    gdf = fix_antimeridian(gdf)
    
    # 2. Convert to Fiji's national metric grid (UTM 58S is good for Fiji)
    # EPSG 32758 is UTM Zone 58S
    gdf = gdf.to_crs(epsg=32758)
    
    # 3. Calculate length_m
    gdf['length_m'] = gdf.geometry.length
    
    # 4. Ensure all are MultiLineString (already handled in fix_antimeridian but just in case)
    new_geoms = []
    for geom in gdf.geometry:
        if geom.geom_type == 'LineString':
            new_geoms.append(MultiLineString([geom]))
        else:
            new_geoms.append(geom)
    gdf.geometry = new_geoms
    
    # 5. Save
    gdf.to_file('fiji_transects_fmg.geojson', driver='GeoJSON')
    print("Successfully saved fiji_transects_fmg.geojson")

if __name__ == "__main__":
    main()