crs-l2-fiji-antimeridian

Rerun

Let me analyze this data carefully.

Looking at the data, I see features that have the `crosses_antimeridian_flag: true`. These are transects whose coordinates cross the 180° meridian (anti-meridian). The issue is that longitudes are stored in the range [-180, 180], but when the transect crosses the anti-meridian, the coordinates jump from ~+179 to ~-179, which would render as an arc spanning the entire globe.

For example, T001 has coordinates:
- [178.94, -17.63]
- [179.46, -17.62]
- [179.97, -17.61]
- [-179.51, -17.60] - this jumps across the antimeridian
- [-178.99, -17.58]
- [-178.48, -17.58]

The fix is to "unwrap" the longitudes for transects that cross the antimeridian - convert the negative longitudes to their positive equivalents (e.g., -179.51 becomes 180.49) so the coordinates flow continuously. Then when we split at 180°, we get two segments that don't span the globe.

Wait, let me think more carefully about the approach.

**Approach:**
1. Read the GeoJSON
2. For each feature:
   - Check if it crosses the antimeridian
   - If it does, "unwrap" the longitudes (shift negative longitudes by +360° so the coordinate sequence is monotonic)
   - Split the unwrapped line at 180° longitude to create a MultiLineString
   - If it doesn't cross the antimeridian, just convert to MultiLineString (with one part)
3. Reproject from WGS84 to Fiji's national metric grid
4. Calculate length in meters
5. Write the output

**Fiji's national metric grid**: Fiji uses two UTM zones (Zone 1 for west of 180°, Zone 60 for east of 180°). Actually, I need to check what the proper Fiji national grid CRS is.

Fiji typically uses:
- **Fiji 1956 / Fiji Map Grid** (EPSG:2136) - a metric grid system for Fiji

Actually, looking up Fiji's national mapping grid:
- EPSG:2136 - Fiji 1956 / Fiji Map Grid (in meters)
- Fiji Map Grid uses the Fiji 1956 datum, with a metric coordinate system

Let me write the script.

For the antimeridian crossing fix:
1. Take the coordinates and check if there's a big jump (> 180° difference between consecutive longitudes)
2. For those that cross, unwrap by adding 360 to negative longitudes
3. Then split the unwrapped line at longitude = 180° (which corresponds to the antimeridian)
4. Create separate LineStrings for each segment
5. Shift the segments that were originally on the other side back by -360° (or just keep them unwrapped for reprojection - actually the CRS conversion should handle that)

Actually, let me think about this differently. The key issue is rendering. I need to:

1. For antimeridian-crossing lines, split them at the antimeridian into multiple segments
2. Each segment stays in its own longitude hemisphere
3. Then convert to MultiLineString

Let me code this up:

For the unwrapping: iterate through coordinates. When we see a jump from ~+180 to ~-180 (or vice versa, difference > 180), that's the antimeridian crossing. We can shift the longitudes on one side by ±360 to make them continuous.

Algorithm:
1. Start with original coordinates
2. Shift coordinates so they're continuous (unwrap): for each pair, if the difference is > 180, add/subtract 360 to subsequent points
3. Find where the unwrapped line crosses 180° longitude
4. Split at those crossing points
5. For each segment, shift back to original range if needed (but keep the coordinates continuous within each segment)

Let me write a Python script for this.

Actually, a simpler approach:
- For each feature that crosses the antimeridian:
  1. Unwrap coordinates (make them continuous by adding 360 to negative longitudes)
  2. Find intersection points with longitude = 180° (or 180° in the unwrapped system)
  3. Split the line at those intersections
  4. The resulting segments will be on two sides: one with longitudes < 180 (original positive longitudes) and one with longitudes > 180 (converted from negative longitudes)
  5. Shift the > 180 segment back by subtracting 360
  6. Combine all segments into a MultiLineString

Let me write this carefully with shapely.

Actually, using shapely's `intersection` and `difference` operations might be cleaner. Let me think...

Better approach using shapely:
1. Create a line from the coordinates
2. Create a "splitter" line along the antimeridian (longitude = 180°)
3. Use shapely's `split` to split the line
4. The result will be geometry parts on either side

Wait, but the coordinates that jump from +179 to -179 create a line that spans the globe through 0°, so splitting at 180° won't help directly with the original coordinates.

Let me reconsider. The issue is that with coordinates in [-180, 180], a line going from [179, -17] to [-179, -17] actually renders as going from 179°E to 179°W going the wrong way (across the Pacific, through the Americas). We want it to render as going from 179°E to 181°E (which is -179°E).

So the approach should be:
1. Unwrap coordinates: if a coordinate's longitude is negative and the previous was positive (and the jump is huge), add 360 to this and all subsequent negative coordinates (or find the right cumulative shift)
2. Now the coordinates are continuous (e.g., 179, 179.5, 180.5 which was -179.5)
3. Split this continuous line at longitude = 180
4. The parts will be: one with longitudes < 180, one with > 180
5. For the > 180 part, subtract 360 to bring back to [-180, 180] range
6. Now all parts are MultiLineString components

Let me code this.

```python
import json
import geopandas as gpd
from shapely.geometry import LineString, MultiLineString, Point, box
from shapely.ops import split, linemerge
import numpy as np
from copy import deepcopy

# Read GeoJSON
with open('/work/fiji_transects_wgs84.geojson', 'r') as f:
    data = json.load(f)

features_out = []

for feature in data['features']:
    props = feature['properties']
    coords = feature['geometry']['coordinates']
    
    # Check if this crosses the antimeridian
    crosses = props.get('crosses_antimeridian_flag', False)
    
    if crosses:
        # Unwrap coordinates: make them continuous
        unwrapped = []
        shift = 0
        prev_lon = coords[0][0]
        for i, (lon, lat) in enumerate(coords):
            if i > 0:
                diff = lon - prev_lon
                if diff > 180:
                    shift -= 360
                elif diff < -180:
                    shift += 360
            unwrapped.append([lon + shift, lat])
            prev_lon = lon
        
        # Create LineString from unwrapped coordinates
        line = LineString(unwrapped)
        
        # Split at antimeridian (180° in original, or 180 + shift in unwrapped)
        # Actually, we split at longitude = 180
        # The antimeridian in the unwrapped system
        # We need to find where the unwrapped line crosses 180, 540, etc.
        
        # Create a vertical line at each crossing
        # Find the range of longitudes
        lons = [c[0] for c in unwrapped]
        min_lon, max_lon = min(lons), max(lons)
        
        # Find which 180-multiples we cross
        # Start from the next 180 multiple above min_lon
        start_m = int(np.ceil(min_lon / 180)) * 180
        end_m = int(np.floor(max_lon / 180)) * 180
        
        split_lons = list(range(int(start_m), int(end_m) + 1, 180))
        
        # Split the line at these longitudes
        parts = [line]
        for split_lon in split_lons:
            splitter = LineString([(split_lon, -90), (split_lon, 90)])
            new_parts = []
            for part in parts:
                try:
                    result = split(part, splitter)
                    new_parts.extend(result.geoms)
                except:
                    new_parts.append(part)
            parts = new_parts
        
        # Now convert parts back to normal longitudes and create lines
        final_parts = []
        for part in parts:
            if part.geom_type == 'LineString' and len(part.coords) >= 2:
                part_coords = []
                for lon, lat in part.coords:
                    # Normalize longitude back to [-180, 180]
                    norm_lon = ((lon + 180) % 360) - 180
                    part_coords.append([norm_lon, lat])
                final_parts.append(LineString(part_coords))
            elif part.geom_type == 'MultiLineString':
                for subpart in part.geoms:
                    if len(subpart.coords) >= 2:
                        part_coords = []
                        for lon, lat in subpart.coords:
                            norm_lon = ((lon + 180) % 360) - 180
                            part_coords.append([norm_lon, lat])
                        final_parts.append(LineString(part_coords))
        
        # Remove degenerate parts
        final_parts = [p for p in final_parts if p.length > 0]
        
        if len(final_parts) == 1:
            geom = final_parts[0]
        else:
            geom = MultiLineString(final_parts)
    else:
        # Just wrap in MultiLineString if not already
        geom = LineString(coords)
    
    features_out.append({
        'type': 'Feature',
        'properties': {
            'transect_id': props['transect_id'],
            'vessel': props['vessel'],
            'survey_date': props['survey_date']
        },
        'geometry': {
            'type': 'MultiLineString' if (crosses and len(final_parts) > 1) or (not crosses) else 'MultiLineString',
            'coordinates': []
        }
    })

# Actually, let me use geopandas instead which handles CRS conversion more cleanly
```

Hmm, this is getting complicated. Let me use a cleaner approach with geopandas and shapely.

Actually, let me use the shapely `split` approach correctly, and then use geopandas for CRS transformation and length calculation.

Let me write the final script:

```python
import json
import geopandas as gpd
from shapely.geometry import LineString, MultiLineString, box
from shapely.ops import split, linemerge
import numpy as np

# Read the GeoJSON
gdf = gpd.read_file('/work/fiji_transects_wgs84.geojson')
gdf = gdf.set_crs('EPSG:4326')

# Fiji Map Grid is EPSG:2136
# But we need to first fix geometries, then reproject

results = []

for idx, row in gdf.iterrows():
    geom = row.geometry
    props = {
        'transect_id': row['transect_id'],
        'vessel': row['vessel'],
        'survey_date': row['survey_date']
    }
    crosses = row.get('crosses_antimeridian_flag', False)
    
    if crosses and geom.geom_type == 'LineString':
        coords = list(geom.coords)
        
        # Unwrap longitudes
        unwrapped = []
        shift = 0
        for i, (lon, lat) in enumerate(coords):
            if i > 0:
                prev_lon = coords[i-1][0]
                diff = lon - prev_lon
                if diff > 180:
                    shift -= 360
                elif diff < -180:
                    shift += 360
            unwrapped.append((lon + shift, lat))
        
        # Create unwrapped line
        unwrapped_line = LineString(unwrapped)
        
        # Find all 180° crossings
        lons = [c[0] for c in unwrapped]
        min_lon, max_lon = min(lons), max(lons)
        
        # Find 180 multiples between min and max
        split_vals = []
        m_start = int(np.ceil(min_lon / 180)) * 180
        m_end = int(np.floor(max_lon / 180)) * 180
        split_vals = list(range(int(m_start), int(m_end) + 1, 180))
        
        # Split at each
        parts = [unwrapped_line]
        for sv in split_vals:
            splitter = LineString([(sv, -90), (sv, 90)])
            new_parts = []
            for part in parts:
                try:
                    result = split(part, splitter)
                    new_parts.extend(result.geoms)
                except:
                    new_parts.append(part)
            parts = new_parts
        
        # Normalize longitudes back
        normalized_parts = []
        for part in parts:
            if part.geom_type == 'LineString' and len(part.coords) >= 2:
                pcoords = []
                for lon, lat in part.coords:
                    norm_lon = ((lon + 180) % 360) - 180
                    pcoords.append((norm_lon, lat))
                normalized_parts.append(LineString(pcoords))
        
        # Filter out degenerate
        normalized_parts = [p for p in normalized_parts if p.length > 1e-10]
        
        if len(normalized_parts) > 1:
            out_geom = MultiLineString(normalized_parts)
        else:
            out_geom = normalized_parts[0]
    else:
        out_geom = geom
    
    # Ensure MultiLineString
    if out_geom.geom_type == 'LineString':
        out_geom = MultiLineString([out_geom])
    
    results.append({'geometry': out_geom, **props})

# Create output GeoDataFrame
out_gdf = gpd.GeoDataFrame(results, geometry='geometry', crs='EPSG:4326')

# Reproject to Fiji Map Grid (EPSG:2136)
out_gdf = out_gdf.to_crs('EPSG:2136')

# Calculate length in meters
out_gdf['length_m'] = out_gdf.geometry.length

# Write to GeoJSON
out_gdf.to_file('/work/fiji_transects_fmg.geojson', driver='GeoJSON')
```

Wait, let me check if EPSG:2136 is the right one. Fiji 1956 / Fiji Map Grid. Let me verify.

Actually, let me check whether EPSG:2136 covers the whole Fiji area, and whether it's the standard national grid for Fiji.

From memory, Fiji Map Grid (FMG) uses the Fiji 1956 datum and is a transverse Mercator projection. Let me use it.

But actually, for some of these lines near 180°, the UTM zones for Fiji are:
- Zone 60 (UTM 60S) covers 174°E to 180°
- Zone 1 (UTM 1S) covers 180° to 186°E (which is -180 to -174°W)

So using a single UTM-based projection like EPSG:2136 (Fiji Map Grid) which is designed for Fiji would be the right approach.

Let me also double-check by looking up EPSG:2136 more carefully... Actually, let me just run the script and see what happens.

I also need to be careful about which parts of the antimeridian-crossing transects get which side. Let me trace through T001:

Original coords: [178.94, -17.63], [179.46, -17.62], [179.97, -17.61], [-179.51, -17.60], [-179.00, -17.58], [-178.48, -17.58]

Unwrapped: 
- [178.94, -17.63] (shift=0)
- [179.46, -17.62] (diff from prev = 0.52, fine)
- [179.97, -17.61] (diff from prev = 0.51, fine)
- [-179.51, -17.60] (diff from prev = -359.48, so shift += 360) → [180.49, -17.60]
- [-179.00, -17.58] (shift=360) → [181.00, -17.58]
- [-178.48, -17.58] (shift=360) → [181.52, -17.58]

So unwrapped: 178.94, 179.46, 179.97, 180.49, 181.00, 181.52

Split at 180: 
- Part 1: 178.94→179.97 (before 180), intersection at 180
  After normalizing: 178.94, 179.46, 179.97 - all fine
- Part 2: 180 → 181.52 (after 180)
  After normalizing: 180 - 360 = -180, 180.49-360 = -179.51, 181.00-360 = -179.00, 181.52-360 = -178.48
  → [-180, -17.60], [-179.51, -17.60], [-179.00, -17.58], [-178.48, -17.58]

But wait, the split should give us an intersection point at exactly 180° longitude. Let me make sure the split operation adds that vertex.

The split function in shapely should add intersection points correctly.

Let me also think about the intersection point normalization. The intersection point is at (180, some_lat). When normalized, 180 → -180 (since ((180+180)%360)-180 = (360%360)-180 = 0-180 = -180). Hmm, that's -180 which is fine.

Actually wait, ((180+180)%360)-180 = 360%360-180 = 0-180 = -180. Okay good.

But 180 is actually on the boundary. Let me think... The point at exactly longitude 180 should stay as 180 or -180. Let me handle this more carefully.

Actually for practical purposes, -180 and 180 are the same line. The important thing is that our parts don't span the globe.

Let me just write and test the script.

Actually, I realize there might be an issue with the split at exactly 180 for the parts. Let me also think about what happens when the intersection point is at 180 and we normalize:

Part 1 intersection: (180, some_lat) → ((180+180)%360)-180 = -180
Part 2 intersection: (180, some_lat) → -180 as well

So both parts share the point at -180 longitude, which connects them correctly.

Let me write the script properly now.

Actually, I just realized I should also handle the case where a transect doesn't cross the antimeridian but the geometry is already a MultiLineString or needs to be converted.

Let me write it all cleanly.

One more thing: the CRS. Let me check what Fiji's national grid CRS is called.

Fiji Map Grid - I believe it's EPSG:2136. Let me just use that.

Let me also think about whether any of the transects that DON'T have the flag set might still have issues. Looking at the data:

T011 through T030 all have `crosses_antimeridian_flag: false` and their coordinates are all on one side of 180°... let me verify a few:
- T011: all longitudes are negative (-176 to -177) - fine, west side
- T013: all negative (-176 to -178) - fine
- T015: all negative (-177 to -178) - fine
- T024: all negative (-179 to -176) - fine but wait, T024 has -179.1 which is very close... let me check if it crosses 180.

T024: [-179.108...], [-178.65...], [-178.19...], [-177.74...], [-177.28...], [-176.82...]
All are between -180 and 0, so no antimeridian crossing. Good.

OK let me write the script.

Actually, hmm, I wonder about the split function. Let me use a different approach: instead of splitting, I can manually cut the unwrapped line at each 180° crossing.

Here's a cleaner approach:

```python
def split_at_antimeridian(coords):
    """Split a line crossing the antimeridian into multiple segments."""
    # Unwrap
    unwrapped = []
    shift = 0
    for i, (lon, lat) in enumerate(coords):
        if i > 0:
            diff = lon - coords[i-1][0]
            if diff > 180:
                shift -= 360
            elif diff < -180:
                shift += 360
        unwrapped.append((lon + shift, lat))
    
    # Find intersection points with multiples of 180
    segments = []
    current_seg = []
    
    for i in range(len(unwrapped) - 1):
        x1, y1 = unwrapped[i]
        x2, y2 = unwrapped[i + 1]
        
        current_seg.append((x1, y1))
        
        # Check if this segment crosses a 180 boundary
        crosses = []
        # Find all 180-multiples between x1 and x2
        start_mlt = int(np.ceil(min(x1, x2) / 180)) * 180
        end_mlt = int(np.floor(max(x1, x2) / 180)) * 180
        
        for mlt in np.arange(start_mlt, end_mlt + 1, 180):
            if min(x1, x2) < mlt < max(x1, x2):
                # Linear interpolation
                t = (mlt - x1) / (x2 - x1)
                y_interp = y1 + t * (y2 - y1)
                crosses.append((mlt, y_interp))
        
        for cp in crosses:
            current_seg.append(cp)
            # Normalize current segment and add
            seg_coords = [((((x+180)%360)-180), y) for x, y in current_seg]
            segments.append(LineString(seg_coords))
            current_seg = [cp]
    
    # Add last point
    current_seg.append(unwrapped[-1])
    seg_coords = [((((x+180)%360)-180), y) for x, y in current_seg]
    segments.append(LineString(seg_coords))
    
    # Filter degenerate
    segments = [s for s in segments if s.length > 1e-10]
    
    return segments
```

This is cleaner and avoids the shapely split function which can be unreliable. Let me finalize the script.

Let me also think about what to do for T010 which has coords: [177.60, -17.28], [179.13, -17.47], [-179.35, -17.65], [-177.82, -17.84]
Unwrapped: 177.60, 179.13, 180.65 (= -179.35+360), 182.18 (= -177.82+360)
Split at 180: segment 1 = 177.60→179.13→180 (interp), segment 2 = 180→180.65→182.18
Norm segment 1: 177.60, 179.13, -180
Norm segment 2: -180, -179.35, -177.82

That looks correct!

OK let me write the final script now.

Let me also handle the case where we might get a 0-length segment or very small segment from the interpolation.

And one more thing: the split_at_antimeridian approach above might not work correctly if there are multiple 180-crossings in one segment of the original line... but looking at the data, each transect with the flag crosses the antimeridian exactly once, so this should be fine.

Actually, I need to be more careful. Let me re-examine the segment accumulation logic.

For T001 unwrapped: 178.94, 179.46, 179.97, 180.49, 181.00, 181.52

Iteration:
i=0: (178.94, 179.46). current_seg = [178.94]. No 180 crossing between 178.94 and 179.46. current_seg = [178.94].
i=1: (179.46, 179.97). current_seg = [178.94, 179.46]. No 180 crossing between 179.46 and 179.97. current_seg = [178.94, 179.46].
i=2: (179.97, 180.49). current_seg = [178.94, 179.46, 179.97]. 180 is between 179.97 and 180.49. 
  Interpolate: t=(180-179.97)/(180.49-179.97)=0.03/0.52=0.0577, y=-17.61+0.0577*(-17.60+17.61)=-17.61+0.000577=-17.6094
  cross point = (180, -17.6094)
  current_seg = [178.94, 179.46, 179.97, 180, -17.6094]
  Norm: [(178.94, -17.63), (179.46, -17.62), (179.97, -17.61), (-180, -17.6094)]
  Segment 1 = LineString([178.94, -17.63], ...)
  current_seg = [180, -17.6094]

i=3: (180.49, 181.00). current_seg = [180, -17.6094, 180.49, -17.60]. No 180 crossing between 180.49 and 181.00 (both > 180). current_seg = [180, -17.6094, 180.49, -17.60]
i=4: (181.00, 181.52). current_seg = [180, -17.6094, 180.49, -17.60, 181.00, -17.58]. No 180 crossing between 181.00 and 181.52. current_seg already contains all remaining points.

After loop: current_seg = [180, -17.6094, 180.49, -17.60, 181.00, -17.58, 181.52, -17.58]
Norm: [-180, -17.6094, -179.51, -17.60, -179.00, -17.58, -178.48, -17.58]
Segment 2 = LineString([-180, -17.6094], [-179.51, -17.60], [-179.00, -17.58], [-178.48, -17.58])

So we get 2 segments:
1. (178.94, -17.63) → (179.46, -17.62) → (179.97, -17.61) → (-180, -17.609)
2. (-180, -17.609) → (-179.51, -17.60) → (-179.00, -17.58) → (-178.48, -17.58)

These form a MultiLineString. When rendered, the first segment goes from 178.94 eastward to -180 (the antimeridian), and the second goes from -180 westward to -178.48. Perfect!

Now let me write the final script.

Let me think about what CRS to use. Let me search for the right EPSG code.

Fiji Map Grid (FMG) - EPSG:2136 (Fiji 1956 / Fiji Map Grid)
This uses the Fiji 1956 datum and transverse Mercator projection.

Actually wait, I should double check that. Let me just try it.

Also, regarding the output: "Carry the survey attributes (vessel and survey_date) through unchanged, use transect_id as the key field"

And: "Every feature must come out as a MultiLineString, so do not leave any transects as plain LineString"

So even non-crossing transects should be wrapped in MultiLineString.

Let me finalize:

```python
import json
import numpy as np
import geopandas as gpd
from shapely.geometry import LineString, MultiLineString, shape

def split_at_antimeridian(coords):
    """
    Split a LineString that crosses the antimeridian (180° longitude)
    into segments that don't cross it.
    Returns a list of LineString segments.
    """
    # Unwrap longitudes to make the coordinate sequence continuous
    unwrapped = []
    shift = 0.0
    for i, (lon, lat) in enumerate(coords):
        if i > 0:
            diff = lon - coords[i-1][0]
            if diff > 180:
                shift -= 360
            elif diff < -180:
                shift += 360
        unwrapped.append((lon + shift, lat))
    
    # Walk through unwrapped segments, splitting at each 180° multiple
    segments = []
    current = []
    
    for i in range(len(unwrapped)):
        x1, y1 = unwrapped[i]
        
        if i == 0:
            current.append((x1, y1))
            continue
        
        x0, y0 = unwrapped[i-1]
        
        # Find all 180 multiples between x0 and x1
        lo, hi = min(x0, x1), max(x0, x1)
        # Multiples of 180 strictly between lo and hi
        start_m = int(np.ceil(lo / 180)) * 180
        end_m = int(np.floor(hi / 180)) * 180
        crossings = list(np.arange(start_m, end_m + 1, 180))
        crossings = [c for c in crossings if lo < c < hi]
        
        if not crossings:
            current.append((x1, y1))
        else:
            # Append points up to each crossing, closing the current segment
            prev_x, prev_y = x0, y0
            for cx in crossings:
                t = (cx - prev_x) / (x1 - prev_x)
                cy = prev_y + t * (y1 - prev_y)
                current.append((cx, cy))
                # Normalize and finalize current segment
                seg_coords = [((((x + 180) % 360) - 180), y) for x, y in current]
                seg = LineString(seg_coords)
                if seg.length > 1e-10:
                    segments.append(seg)
                # Start new segment from this crossing point
                current = [(cx, cy)]
                prev_x, prev_y = cx, cy
            
            # Add the current endpoint
            current.append((x1, y1))
    
    # Finalize last segment
    if len(current) >= 2:
        seg_coords = [((((x + 180) % 360) - 180), y) for x, y in current]
        seg = LineString(seg_coords)
        if seg.length > 1e-10:
            segments.append(seg)
    
    return segments


# Read input
with open('/work/fiji_transects_wgs84.geojson') as f:
    data = json.load(f)

features_out = []

for feat in data['features']:
    props = feat['properties']
    geom = feat['geometry']
    coords = geom['coordinates']
    crosses = props.get('crosses_antimeridian_flag', False)
    
    if crosses and geom['type'] == 'LineString':
        parts = split_at_antimeridian(coords)
        if len(parts) == 0:
            continue  # skip degenerate
        out_geom = MultiLineString(parts) if len(parts) > 1 else parts[0]
    else:
        out_geom = shape(geom)
    
    # Ensure MultiLineString
    if out_geom.geom_type == 'LineString':
        out_geom = MultiLineString([out_geom])
    
    features_out.append({
        'type': 'Feature',
        'properties': {
            'transect_id': props['transect_id'],
            'vessel': props['vessel'],
            'survey_date': props['survey_date']
        },
        'geometry': {
            'type': 'MultiLineString',
            'coordinates': [
                [list(c) for c in part.coords]
                for part in out_geom.geoms
            ]
        }
    })

# Write intermediate fixed GeoJSON
fc = {
    'type': 'FeatureCollection',
    'features': features_out
}
with open('/work/fiji_fixed_wgs84.geojson', 'w') as f:
    json.dump(fc, f)

# Now use geopandas to reproject and add length
gdf = gpd.read_file('/work/fiji_fixed_wgs84.geojson')
gdf = gdf.set_crs('EPSG:4326')
gdf = gdf.to_crs('EPSG:2136')
gdf['length_m'] = gdf.geometry.length.round(2)
gdf.to_file('/work/fiji_transects_fmg.geojson', driver='GeoJSON')
```

OK, I think this approach is solid. Let me code it up properly and test.

Wait, I realized I should just use geopandas throughout instead of manually creating the output geojson. Let me simplify.

Actually, the manual approach for the splitting is needed since geopandas doesn't have a built-in antimeridian splitter. Let me keep the split function and then use geopandas for the CRS conversion and length calculation.

Let me finalize and write the script.

Good, the script ran without errors. Let me verify the output now.

1. Check it's valid GeoJSON
2. Check all features are MultiLineString
3. Check transect_id, vessel, survey_date are preserved
4. Check length_m is present
5. Check the CRS is EPSG:2136 (Fiji Map Grid) 
6. Check the antimeridian-crossing transects are properly split

Wait, something is wrong. The lat values are around -60,000,000 to -61,000,000 meters. That's way too far south for Fiji. Let me check what EPSG:2136 actually is.

EPSG:2136 is "Fiji 1956 / Fiji Map Grid" - a projected CRS for Fiji. Let me check the standard parameters.

Actually, looking at the coordinates, those are huge negative numbers in the millions. That's not right for a metric grid. Fiji is in the southern hemisphere so y (northing) should be positive if using a standard false northing, or negative if false northing is 0. But -60 million meters is extreme.

Let me check - EPSG:2136 uses the "Fiji 1956" datum, which is a local datum. The projection is likely a Transverse Mercator with a central meridian of 177°E or 178°E or 180°... Let me look more carefully.

Actually, EPSG:2136 - Fiji 1956 / Fiji Map Grid. The projection is "Fiji Map Grid" which is a Transverse Mercator projection with:
- Central Meridian = 178°E (some sources say 180°...)

Let me check by looking at the actual EPSG definition. Actually, I think the issue is that EPSG:2136 uses a Hotine Oblique Mercator projection, not a standard Transverse Mercator. The definition is quite specific.

Wait, let me reconsider. The lat values around -60 million suggest something is off. For a location around 18°S in Fiji, the y coordinate should be something like:
- If false northing is 0, then y ≈ -18° in meters = -18 * 111,000 ≈ -2,000,000 meters
- But -60,000,000 is way off.

Actually, maybe I'm using the wrong EPSG code. Let me think about what "Fiji's national metric grid" actually is.

The modern Fiji national grid is EPSG:2136 "Fiji 1956 / Fiji Map Grid", but there's also:
- EPSG:3460 "Fiji 1986 / Fiji Map Grid"
- Or maybe they use a more recent CRS

Wait, let me look at the longitude values more carefully. For T011 which is at longitudes around -176.8 to -177.6 degrees, the projected x values are -539,630 to -284,109. That seems plausible for an Oblique Mercator projection centered on Fiji.

But the y values at -60,800,000 to -61,200,000 are way off. These look like they might be in millimeters rather than meters, or something else is wrong.

Actually, I think the issue might be that pyproj / geopandas is using a different datum transformation and the coordinates are computed incorrectly. Let me check what pyproj says about EPSG:2136.

Actually, let me look at this more carefully. The Fiji Map Grid coordinates should look like:
- Easting: around 100,000 to 900,000 meters (typical for a projected CRS with a false easting)
- Northing: positive values in the southern hemisphere, or negative if false northing is 0

For Oblique Mercator projections, the coordinates can be quite different from UTM. Let me check.

Actually, I realize the y values are about -60 million, which is about the right magnitude if you consider degrees * 111,111 * 3 = degrees * 333,333... No, that doesn't make sense.

Let me actually check the EPSG definition parameters for 2136.

From the EPSG database, Fiji 1956 / Fiji Map Grid (EPSG:2136) is:
- Method: Hotine Oblique Mercator (variant A)
- Latitude of projection center: -17° 50' (approx 17.8333°S)
- Longitude of projection center: 178° 20' (178.3333°E) 
- Azimuth: 332.5°
- Scale factor: 1.0
- False easting: 300,000
- False northing: 500,000

Hmm, but the coordinates I'm getting look very wrong. Let me just run a quick check with pyproj to see what's happening.

Actually, let me just try using EPSG:3460 (Fiji 1986 / Fiji Map Grid) instead, or maybe I should look at what the actual coordinates should be.

Wait, let me reconsider. The pyproj transformation from WGS84 to Fiji 1956 might be using a bad transformation or no transformation. Let me check by reprojecting a single point.

Ah, the issue is clear. EPSG:2136 is not Fiji Map Grid - it's Accra / Ghana National Grid. Fiji Map Grid is EPSG:3460. Let me fix the CRS.

Also, looking at the coordinate output from the current file, the lat values are around -60,000,000 which is clearly wrong for Fiji (which should be around 3,000,000 to 4,000,000 northing for FMG). The eastings should be around 2,000,000. So clearly EPSG:2136 was wrong.

Let me also look more carefully at the antimeridian splitting. The coordinates with longitudes crossing from positive (east of 180) to negative (west of 180) represent transects that cross the 180° meridian (antimeridian). The issue is that in WGS84, coordinates like [179.97..., -17.61...] then [-179.51..., -17.60...] appear to jump across the whole globe (from near 180°E to near 180°W), when really they're just crossing the antimeridian.

For the Fiji Map Grid (EPSG:3460), the central meridian is 178.75°E, so this should handle Fiji data well.

Let me fix the script to use EPSG:3460 instead of EPSG:2136.