crs-l2-fiji-antimeridian

Rerun

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

def fix_antimeridian_lines(gdf):
    """
    Fixes LineStrings that cross the antimeridian by splitting them into 
    multiple segments and wrapping longitudes to a consistent range.
    """
    new_features = []
    
    for _, row in gdf.iterrows():
        geom = row.geometry
        if geom is None:
            continue
            
        # Get coordinates
        coords = list(geom.coords)
        
        # Check for antimeridian crossing (jump in longitude > 180)
        # However, in WGS84, the jump is usually from +180 to -180 or vice versa.
        # A simple way is to detect if the diff between consecutive longitudes is > 180.
        
        split_points = []
        for i in range(len(coords) - 1):
            lon1, lat1 = coords[i]
            lon2, lat2 = coords[i+1]
            
            diff = lon2 - lon1
            if diff > 180:
                # Jump from -179 to 179 (wrong direction for a short segment)
                # or jump from 179 to -179 (correct direction)
                # Actually, if diff > 180, it means it's going "the long way" 
                # around the world if we assume it should be a short segment.
                # But wait, if lon1=179 and lon2=-179, diff = -358.
                # If lon1=-179 and lon2=179, diff = 358.
                pass
                
            # Let's use a more robust logic:
            # If the absolute difference is > 180, it crossed the antimeridian.
            if abs(diff) > 180:
                # We need to split at this point.
                # But where is the "split" point? The segment itself crosses it.
                # For a LineString, we can't easily "split" a single segment unless we 
                # interpolate the crossing point. 
                # However, the prompt says "render as arcs spanning the entire globe".
                # This usually means the coordinates are like [179, -17], [-179, -17].
                # When plotted, this goes the long way around.
                # If we want it to be a short segment, we should treat the jump 
                # as crossing the antimeridian.
                
                # A better approach:
                # 1. Normalize all longitudes to be continuous.
                # 2. Detect where the jump occurs.
                # 3. Split the geometry there.
                pass

        # Let's try a different approach:
        # For each segment, if abs(lon2 - lon1) > 180:
        #   adjust lon2 so that |lon2 - lon1| < 180
        #   this will create a "jump" in the coordinates, which we can then split.
        
        adjusted_coords = []
        if len(coords) > 0:
            adjusted_coords.append(coords[0])
            for i in range(len(coords) - 1):
                lon1, lat1 = coords[i]
                lon2, lat2 = coords[i+1]
                
                diff = lon2 - lon1
                if diff > 180:
                    lon2 -= 360
                elif diff < -180:
                    lon2 += 360
                
                adjusted_coords.append((lon2, lat2))
        
        # Now adjusted_coords might have values like 181, 182...
        # We want to split the LineString where the longitude crosses the 180/-180 boundary.
        # Or more simply, where the "wrapped" longitude would cross 180 or -180.
        
        # Let's create a LineString from adjusted_coords.
        # Then we split it where the longitude is 180 or -180.
        # But we need to be careful with the coordinates.
        
        # Actually, the most robust way to handle antimeridian in GeoPandas/Shapely
        # is to use the 'antimeridian' library if available, but it's not.
        # Let's do it manually:
        
        if not adjusted_coords:
            continue

        # Create segments
        segments = []
        current_segment = [adjusted_coords[0]]
        
        for i in range(1, len(adjusted_coords)):
            prev_lon = adjusted_coords[i-1][0]
            curr_lon = adjusted_coords[i][0]
            
            # Check if this segment crosses the 180 line
            # We normalize the line to a continuous range, then split at 180.
            # But the "jump" happens when we transition between the 180/-180 boundary.
            
            # Let's check if the segment crosses the boundary (wrapping around 180)
            # Since we adjusted lon2, a crossing happens if the segment 
            # crosses an integer multiple of 180.
            
            # Actually, let's just use the adjusted_coords and split 
            # when the coordinate crosses the boundary.
            
            # Wait, if I have coords [179, 181], it's a segment crossing the line.
            # I can split it at 180.
            
            # To keep it simple:
            # 1. Generate segments that don't cross the 180 meridian.
            # 2. For each segment, normalize its coordinates back to [-180, 180].
            
            # Let's refine the "split" logic.
            # For each segment (p1, p2) in adjusted_coords:
            #   if the segment crosses a multiple of 180:
            #     split it into (p1, p_boundary) and (p_boundary, p2)
            
            # But wait, if we just want to FIX the "arcs", we just need to 
            # make sure the line doesn't go the "long way".
            # The "long way" is when |lon2 - lon1| > 180.
            # The "short way" is when |lon2 - lon1| < 180.
            # My `adjusted_coords` logic already converts "long way" to "short way"
            # by adjusting lon2.
            # For example: [179, -179] -> diff = -358. diff < -180, so lon2 += 360 -> lon2 = 181.
            # Coords: [179, -179] -> [179, 181]. This is a short segment.
            # Now we just need to split this LineString where it crosses 180.
            
            # Let's do that.
            pass
            
        # Re-implementing the logic inside the loop properly.
        pass

def solve():
    gdf = gpd.read_file('fiji_transects_wgs84.geojson')
    
    # 1. Fix Antimeridian
    fixed_geoms = []
    for _, row in gdf.iterrows():
        geom = row.geometry
        if geom is None or geom.is_empty:
            fixed_geoms.append(geom)
            continue
            
        coords = list(geom.coords)
        adjusted_coords = [coords[0]]
        for i in range(len(coords) - 1):
            lon1, lat1 = coords[i]
            lon2, lat2 = coords[i+1]
            diff = lon2 - lon1
            if diff > 180:
                lon2 -= 360
            elif diff < -180:
                lon2 += 360
            adjusted_coords.append((lon2, lat2))
            
        # Now we have a continuous path in "unwrapped" longitude space.
        # We need to split this path whenever it crosses a multiple of 180.
        # Or more accurately, whenever the longitude crosses the boundary.
        # Since we want the output in [-180, 180], we split at 180, -180, etc.
        
        # A simpler way to split:
        # Iterate through segments. If a segment crosses a boundary, 
        # split it and start a new segment.
        
        sub_segments = []
        current_sub_coords = [adjusted_coords[0]]
        
        for i in range(1, len(adjusted_coords)):
            p1 = adjusted_coords[i-1]
            p2 = adjusted_coords[i]
            
            # Check if the segment [p1, p2] crosses a multiple of 180
            # The multiples of 180 are ... -180, 0, 180, 360 ...
            # But we only care about the boundaries of the [-180, 180] range.
            # Since we are in "unwrapped" space, we look for crossings of 180, -180, 360, etc.
            
            # Wait, the simplest way to find if it crosses is to check if 
            # floor(lon1/180) != floor(lon2/180) ? No, that's not quite right.
            # The boundaries are ..., -180, 180, 360, ...
            # If the segment crosses 180, or -180, or 360...
            
            # Let's find if there is an integer k such that 180*k is between lon1 and lon2.
            # (Wait, 0 is also a multiple, but we don't want to split at 0).
            # The boundaries are 180, -180, 360, -360, etc.
            # So we check for crossings of 180 * k where k is an ODD integer.
            # 180 * 1 = 180
            # 180 * -1 = -180
            # 180 * 3 = 540
            # 180 * -3 = -540
            
            # Actually, even simpler: 
            # In our "unwrapped" space, the boundaries are at ... -180, 180, 360 ...
            # These are (2k+1)*180.
            
            # Let's check if the segment (p1, p2) crosses any (2k+1)*180.
            
            def get_boundary(lon1, lon2):
                # Find k such that (2k+1)*180 is between lon1 and lon2
                # We can check if floor(lon/180) is odd or even? 
                # Let's just check the range.
                
                # Check for the crossing point
                # The boundaries are 180, -180, 360, -360, 540, -540...
                # We can iterate through possible boundaries.
                # But since the segments are short, there's at most one.
                
                # Let's try to find the boundary.
                # We are looking for B = (2k+1)*180
                # such that min(lon1, lon2) < B < max(lon1, lon2)
                
                # We can estimate k.
                # B/180 = 2k+1 => k = (B/180 - 1) / 2
                # Let's just check a few k's around lon1/180.
                
                low_k = int(lon1 / 180) - 2
                high_k = int(lon1 / 180) + 2
                for k in range(low_k, high_k + 1):
                    B = (2 * k + 1) * 180
                    if min(lon1, lon2) < B < max(lon1, lon2):
                        # Found it!
                        # Interpolate to find the exact lon at this B
                        # (Actually, for a LineString, we usually just split at the boundary)
                        # We need the lat at this B.
                        # Linear interpolation:
                        t = (B - lon1) / (lon2 - lon1)
                        lat_at_B = lat1 + t * (lat2 - lat1)
                        return B, lat_at_B
                return None

            boundary = get_boundary(p1[0], p2[0])
            if boundary:
                B_lon, B_lat = boundary
                # Split: current_sub_coords ends at (B_lon, B_lat)
                # New segment starts at (B_lon, B_lat)
                current_sub_coords.append((B_lon, B_lat))
                sub_segments.append(LineString(current_sub_coords))
                current_sub_coords = [(p2[0], p2[1])] # This is wrong, should be the second part of the split
                # Wait, the split is: segment is (p1, p2). 
                # We split it at B into (p1, B) and (B, p2).
                # But our `current_sub_coords` already has p1.
                # So the segment is [..., p1, B].
                # The next segment is [B, p2].
                # However, we are iterating segment by segment.
                # If we find a boundary in (p1, p2), we:
                # 1. finish current_sub_coords with (B_lon, B_lat).
                # 2. start a new current_sub_coords with (B_lon, B_lat) and (p2[0], p2[1]).
                
                # Let's redo this loop logic.
                pass
        
        # Corrected loop logic:
        sub_segments = []
        current_sub_coords = [adjusted_coords[0]]
        for i in range(1, len(adjusted_coords)):
            p1 = adjusted_coords[i-1]
            p2 = adjusted_coords[i]
            
            boundary = get_boundary(p1[0], p2[0])
            if boundary:
                B_lon, B_lat = boundary
                # Finish current sub-segment
                current_sub_coords.append((B_lon, B_lat))
                sub_segments.append(LineString(current_sub_coords))
                # Start new sub-segment
                current_sub_coords = [(B_lon, B_lat), (p2[0], p2[1])]
            else:
                current_sub_coords.append(p2)
        
        # Last one
        if len(current_sub_coords) >= 2:
            sub_segments.append(LineString(current_sub_coords))
            
        # Now, normalize all coordinates in all sub_segments back to [-180, 180]
        final_sub_geoms = []
        for sg in sub_segments:
            # Normalize longitudes: (lon + 180) % 360 - 180
            new_coords = []
            for x, y in sg.coords:
                new_x = (x + 180) % 360 - 180
                # Handle the case where % 360 returns 0 for -180 (which is fine)
                # but we want to ensure we don't get +180 if we wanted -180. 
                # Actually, -180 and 180 are the same.
                new_coords.append((new_x, y))
            final_sub_geoms.append(LineString(new_coords))
            
        if not final_sub_geoms:
            fixed_geoms.append(geom) # Fallback
        else:
            fixed_geoms.append(MultiLineString(final_sub_geoms))

    gdf.geometry = fixed_geoms

    # 2. Convert to Fiji National Metric Grid
    # Fiji uses Fiji Grid (EPSG:2978 or similar? Let's check)
    # Actually, common for Fiji is Fiji Transverse Mercator.
    # Let's search for Fiji CRS.
    # EPSG:2978 is "Fiji Grid".
    # Let's try to use EPSG:2978.
    
    fiji_crs = "EPSG:2978"
    gdf = gdf.to_crs(fiji_crs)
    
    # 3. Calculate length in metres
    gdf['length_m'] = gdf.geometry.length
    
    # 4. Ensure all are MultiLineString
    def ensure_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom
        
    gdf.geometry = gdf.geometry.map(ensure_multilinestring)
    
    # 5. Save to GeoJSON
    # The prompt asks for 'fiji_transects_fmg.geojson'
    gdf.to_file('fiji_transects_fmg.geojson', driver='GeoJSON')

# ... wait, I need to define get_boundary inside solve or pass it.
# And I'll use a more robust way to ensure MultiLineString.
solve()
` , file_path: `/work/solve.py`

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

def get_boundary(lon1, lon2):
    """
    Find k such that (2k+1)*180 is between lon1 and lon2.
    Returns (B_lon, B_lat) or None.
    """
    # Boundaries are ..., -180, 180, 360, ...
    # We look for B = (2k+1)*180
    
    # We can estimate k.
    # B/180 = 2k+1 => k = (B/180 - 1) / 2
    # Let's check a few k's around lon1/180.
    
    low_k = int(lon1 / 180) - 2
    high_k = int(lon1 / 180) + 2
    for k in range(low_k, high_k + 1):
        B = (2 * k + 1) * 180
        if min(lon1, lon2) < B < max(lon1, lon2):
            return B
    return None

def solve():
    gdf = gpd.read_file('fiji_transects_wgs84.geojson')
    
    # 1. Fix Antimeridian
    fixed_geoms = []
    for _, row in gdf.iterrows():
        geom = row.geometry
        if geom is None or geom.is_empty:
            fixed_geoms.append(geom)
            continue
            
        coords = list(geom.coords)
        adjusted_coords = [coords[0]]
        for i in range(len(coords) - 1):
            lon1, lat1 = coords[i]
            lon2, lat2 = coords[i+1]
            diff = lon2 - lon1
            if diff > 180:
                lon2 -= 360
            elif diff < -180:
                lon2 += 360
            adjusted_coords.append((lon2, lat2))
            
        sub_segments = []
        current_sub_coords = [adjusted_coords[0]]
        
        for i in range(1, len(adjusted_coords)):
            p1 = adjusted_coords[i-1]
            p2 = adjusted_coords[i]
            
            B_lon = get_boundary(p1[0], p2[0])
            if B_lon is not None:
                # Interpolate lat at B_lon
                t = (B_lon - p1[0]) / (p2[0] - p1[0])
                B_lat = p1[1] + t * (p2[1] - p1[1])
                
                # Finish current sub-segment
                current_sub_coords.append((B_lon, B_lat))
                sub_segments.append(LineString(current_sub_coords))
                # Start new sub-segment
                current_sub_coords = [(B_lon, B_lat), (p2[0], p2[1])]
            else:
                current_sub_coords.append(p2)
        
        if len(current_sub_coords) >= 2:
            sub_segments.append(LineString(current_sub_coords))
            
        # Normalize all coordinates in all sub_segments back to [-180, 180]
        final_sub_geoms = []
        for sg in sub_segments:
            new_coords = []
            for x, y in sg.coords:
                # Normalize longitude to [-180, 180]
                # Using (x + 180) % 360 - 180 works well.
                nx = (x + 180) % 360 - 180
                # If nx is -180, it's fine.
                new_coords.append((nx, y))
            final_sub_geoms.append(LineString(new_coords))
            
        if not final_sub_geoms:
            fixed_geoms.append(geom)
        else:
            fixed_geoms.append(MultiLineString(final_sub_geoms))

    gdf.geometry = fixed_geoms

    # 2. Convert to Fiji National Metric Grid (EPSG:2978)
    # Note: EPSG:2978 is Fiji Grid.
    fiji_crs = "EPSG:2978"
    gdf = gdf.to_crs(fiji_crs)
    
    # 3. Calculate length in metres
    gdf['length_m'] = gdf.geometry.length
    
    # 4. Ensure all are MultiLineString
    def ensure_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom
        
    gdf.geometry = gdf.geometry.map(ensure_multilinestring)
    
    # 5. Save to GeoJSON
    # Keeping original columns and adding length_m
    # Columns to keep: transect_id, vessel, survey_date, length_m
    # We'll keep all and ensure the ones requested are there.
    
    # The user said: "Carry the survey attributes (vessel and survey_date) through unchanged, 
    # use transect_id as the key field".
    
    gdf.to_file('fiji_transects_fmg.geojson', driver='GeoJSON')

if __name__ == "__main__":
    solve()

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

def get_boundary(lon1, lon2):
    """
    Find k such that (2k+1)*180 is between lon1 and lon2.
    Returns B_lon or None.
    """
    low_k = int(lon1 / 180) - 2
    high_k = int(lon1 / 180) + 2
    for k in range(low_k, high_k + 1):
        B = (2 * k + 1) * 180
        if min(lon1, lon2) < B < max(lon1, lon2):
            return B
    return None

def solve():
    input_file = 'fiji_transects_wgs84.geojson'
    output_file = 'fiji_transects_fmg.geojson'
    
    if not os.path.exists(input_file):
        print(f"Error: {input_file} not found")
        return

    gdf = gpd.read_file(input_file)
    
    # 1. Fix Antimeridian
    fixed_geoms = []
    for _, row in gdf.iterrows():
        geom = row.geometry
        if geom is None or geom.is_empty:
            fixed_geoms.append(geom)
            continue
            
        coords = list(geom.coords)
        adjusted_coords = [coords[0]]
        for i in range(len(coords) - 1):
            lon1, lat1 = coords[i]
            lon2, lat2 = coords[i+1]
            diff = lon2 - lon1
            if diff > 180:
                lon2 -= 360
            elif diff < -180:
                lon2 += 360
            adjusted_coords.append((lon2, lat2))
            
        sub_segments = []
        current_sub_coords = [adjusted_coords[0]]
        
        for i in range(1, len(adjusted_coords)):
            p1 = adjusted_coords[i-1]
            p2 = adjusted_coords[i]
            
            B_lon = get_boundary(p1[0], p2[0])
            if B_lon is not None:
                t = (B_lon - p1[0]) / (p2[0] - p1[0])
                B_lat = p1[1] + t * (p2[1] - p1[1])
                
                current_sub_coords.append((B_lon, B_lat))
                sub_segments.append(LineString(current_sub_coords))
                current_sub_coords = [(B_lon, B_lat), (p2[0], p2[1])]
            else:
                current_sub_coords.append(p2)
        
        if len(current_sub_coords) >= 2:
            sub_segments.append(LineString(current_sub_coords))
            
        final_sub_geoms = []
        for sg in sub_segments:
            new_coords = []
            for x, y in sg.coords:
                nx = (x + 180) % 360 - 180
                new_coords.append((nx, y))
            final_sub_geoms.append(LineString(new_coords))
            
        if not final_sub_geoms:
            fixed_geoms.append(geom)
        else:
            fixed_geoms.append(MultiLineString(final_sub_geoms))

    gdf.geometry = fixed_geoms

    # 2. Convert to Fiji National Metric Grid (EPSG:2978)
    fiji_crs = "EPSG:2978"
    gdf = gdf.to_crs(fiji_crs)
    
    # 3. Calculate length in metres
    gdf['length_m'] = gdf.geometry.length
    
    # 4. Ensure all are MultiLineString
    def ensure_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom
        
    gdf.geometry = gdf.geometry.map(ensure_multilinestring)
    
    # 5. Save to GeoJSON
    gdf.to_file(output_file, driver='GeoJSON')
    print(f"Successfully saved to {output_file}")

if __name__ == "__main__":
    solve()

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

def get_boundary(lon1, lon2):
    low_k = int(lon1 / 180) - 2
    high_k = int(lon1 / 180) + 2
    for k in range(low_k, high_k + 1):
        B = (2 * k + 1) * 180
        if min(lon1, lon2) < B < max(lon1, lon2):
            return B
    return None

def solve():
    input_file = 'fiji_transects_wgs84.geojson'
    output_file = 'fiji_transects_fmg.geojson'
    
    if not os.path.exists(input_file):
        print(f"Error: {input_file} not found")
        return

    gdf = gpd.read_file(input_file)
    
    fixed_geoms = []
    for _, row in gdf.iterrows():
        geom = row.geometry
        if geom is None or geom.is_empty:
            fixed_geoms.append(geom)
            continue
            
        coords = list(geom.coords)
        adjusted_coords = [coords[0]]
        for i in range(len(coords) - 1):
            lon1, lat1 = coords[i]
            lon2, lat2 = coords[i+1]
            diff = lon2 - lon1
            if diff > 180:
                lon2 -= 360
            elif diff < -180:
                lon2 += 360
            adjusted_coords.append((lon2, lat2))
            
        sub_segments = []
        current_sub_coords = [adjusted_coords[0]]
        
        for i in range(1, len(adjusted_coords)):
            p1 = adjusted_coords[i-1]
            p2 = adjusted_coords[i]
            
            B_lon = get_boundary(p1[0], p2[0])
            if B_lon is not None:
                t = (B_lon - p1[0]) / (p2[0] - p1[0])
                B_lat = p1[1] + t * (p2[1] - p1[1])
                current_sub_coords.append((B_lon, B_lat))
                sub_segments.append(LineString(current_sub_coords))
                current_sub_coords = [(B_lon, B_lat), (p2[0], p2[1])]
            else:
                current_sub_coords.append(p2)
        
        if len(current_sub_coords) >= 2:
            sub_segments.append(LineString(current_sub_coords))
            
        final_sub_geoms = []
        for sg in sub_segments:
            new_coords = []
            for x, y in sg.coords:
                nx = (x + 180) % 360 - 180
                new_coords.append((nx, y))
            final_sub_geoms.append(LineString(new_coords))
            
        if not final_sub_geoms:
            fixed_geoms.append(geom)
        else:
            fixed_geoms.append(MultiLineString(final_sub_geoms))

    gdf.geometry = fixed_geoms
    fiji_crs = "EPSG:2978"
    gdf = gdf.to_crs(fiji_crs)
    gdf['length_m'] = gdf.geometry.length
    
    def ensure_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom
        
    gdf.geometry = gdf.geometry.map(ensure_multilinestring)
    gdf.to_file(output_file, driver='GeoJSON')
    print(f"Successfully saved to {output_file}")

if __name__ == "__main__":
    solve()