crs-l2-fiji-antimeridian

Rerun

import geopandas as gpd

def inspect_geojson(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("\nCRS:", gdf.crs)
    print("\nGeometry types:", gdf.geometry.type.unique())
    print("\nFirst few rows:")
    print(gdf.head())
    print("\nGeometry sample (to check for long arcs):")
    for idx, row in gdf.head(5).iterrows():
        print(f"ID: {row.get('transect_id', idx)}, Type: {row.geometry.type}, Bounds: {row.geometry.bounds}")

if __name__ == "__main__":
    inspect_geojson('fiji_transects_wgs84.geojson')

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

def fix_antimeridian_crossing(line, lon_boundary=180.0):
    """
    Splits a LineString that crosses the antimeridian into multiple LineStrings.
    """
    coords = list(line.coords)
    new_segments = []
    current_segment = []

    for i in range(len(coords)):
        if i > 0:
            # Check for jump across antimeridian
            # A jump is considered if the difference between consecutive longitudes
            # is very large (e.g., > 180)
            diff = coords[i][0] - coords[i-1][0]
            
            if abs(diff) > 180:
                # Close current segment
                if current_segment:
                    new_segments.append(LineString(current_segment))
                    current_segment = []
                
                # Handle the jump
                # If diff > 0, it went from -180 area to +180 area
                # If diff < 0, it went from +180 area to -180 area
                # We adjust the longitude to be continuous in a local sense
                if diff > 0:
                    # Jumped from e.g. 179 to -179 (actually the diff would be -358)
                    # Wait, if coords[i] is -179 and coords[i-1] is 179, diff is -358
                    # If diff < -180, it's a jump from positive to negative (wrapping around)
                    # We add 360 to the new coord to make it continuous for calculation, 
                    # but we'll actually split it.
                    pass 

        # Actually, a better way for GeoJSON is to detect the jump and split.
        # If we detect a jump, we start a new segment.
        # To handle the wrapping, we can normalize the longitudes relative to the first point.
        pass

    # Let's rethink.
    pass

def split_antimeridian_geometry(geom, lon_boundary=180.0):
    """
    Correctly splits a LineString crossing the antimeridian.
    """
    if geom.geom_type != 'LineString':
        return geom

    coords = list(geom.coords)
    segments = []
    current_coords = [coords[0]]

    for i in range(1, len(coords)):
        prev_lon = coords[i-1][0]
        curr_lon = coords[i][0]
        
        # Check if it crosses the antimeridian
        # The jump is detected if the distance between longitudes is > 180
        if abs(curr_lon - prev_lon) > 180:
            # Split here
            segments.append(LineString(current_coords))
            current_coords = [coords[i]]
            
            # Normalize the new segment's longitude relative to the previous one
            # to keep it "close" in a coordinate sense, though for the final 
            # GeoJSON we want the actual WGS84 lon.
            # However, a LineString that crosses the antimeridian is problematic.
            # The task says "fix the geometries so they render correctly".
            # Rendering correctly in most GIS software means splitting the line
            # at the antimeridian.
        else:
            current_coords.append(coords[i])
    
    if current_coords:
        segments.append(LineString(current_coords))
        
    if len(segments) == 1:
        return segments[0]
    else:
        return MultiLineString(segments)

def process_transects(input_file, output_file):
    gdf = gpd.read_file(input_file)
    
    # 1. Fix antimeridian crossings
    # We apply a function that returns a LineString or MultiLineString
    gdf['geometry'] = gdf['geometry'].apply(lambda x: split_antimeridian_geometry(x))
    
    # 2. Convert to Fiji's national metric grid (FMG)
    # Searching for Fiji projected CRS. EPSG:2978 is Fiji Red Beach (often used).
    # Let's try to find a suitable one or use a common one. 
    # Actually, the user says "Fiji's national metric grid". 
    # Usually, this refers to Fiji Grid (EPSG:2978 or similar). 
    # I'll use EPSG:2978 for Fiji Red Beach if available, or a standard UTM zone for Fiji.
    # Fiji is roughly between 160E and 180E.
    # Let's use EPSG:2978 (Fiji Red Beach) if it's the standard.
    # If not, I'll use a generic projection or check.
    # Actually, EPSG:2978 is indeed Fiji Red Beach.
    target_crs = "EPSG:2978"
    
    gdf = gdf.to_crs(target_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'].apply(ensure_multilinestring)
    
    # 5. Save to GeoJSON
    # Keep only required columns (plus transect_id, vessel, survey_date, length_m)
    # The prompt says "Carry the survey attributes (vessel and survey_date) through unchanged"
    # and "use transect_id as the key field".
    # It doesn't explicitly say to discard 'crosses_antimeridian_flag'. 
    # Usually, it's safer to keep all existing attributes unless told otherwise.
    
    gdf.to_file(output_file, driver='GeoJSON')
    print(f"Saved to {output_file}")

if __name__ == "__main__":
    process_transects('fiji_transects_wgs84.geojson', 'fiji_transects_fmg.geojson')

import geopandas as gpd

def inspect_output(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("\nCRS:", gdf.crs)
    print("\nGeometry types:", gdf.geometry.type.unique())
    print("\nFirst few rows:")
    print(gdf.head())
    print("\nChecking for length_m attribute:")
    print(gdf[['transect_id', 'length_m']].head())
    print("\nChecking geometry types specifically:")
    print(gdf.geometry.apply(lambda x: x.geom_type).value_counts())

if __name__ == "__main__":
    inspect_output('fiji_transects_fmg.geojson')

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

def split_antimeridian_geometry(geom):
    """
    Correctly splits a LineString crossing the antimeridian.
    The input geom is in WGS84 (EPSG:4326).
    """
    if geom.geom_type != 'LineString':
        return geom

    coords = list(geom.coords)
    segments = []
    current_coords = [coords[0]]

    for i in range(1, len(coords)):
        prev_lon = coords[i-1][0]
        curr_lon = coords[i][0]
        
        # Check if it crosses the antimeridian
        # The jump is detected if the absolute difference is > 180
        if abs(curr_lon - prev_lon) > 180:
            # Split here
            if current_coords:
                segments.append(LineString(current_coords))
            current_coords = [coords[i]]
        else:
            current_coords.append(coords[i])
    
    if current_coords:
        segments.append(LineString(current_coords))
        
    if len(segments) == 1:
        return segments[0]
    else:
        return MultiLineString(segments)

def ensure_multilinestring(geom):
    """
    Ensures all geometries are MultiLineString.
    """
    if geom.geom_type == 'LineString':
        return MultiLineString([geom])
    elif geom.geom_type == 'MultiLineString':
        return geom
    else:
        # If it's something else, we'll try to handle it, but based on input it's LineString
        return geom

def process_transects(input_file, output_file):
    # Read input
    gdf = gpd.read_file(input_file)
    
    # 1. Fix antimeridian crossings in WGS84
    # We must do this BEFORE projecting to metric, because the projection 
    # will interpret the long arcs as real paths.
    gdf['geometry'] = gdf['geometry'].apply(split_antimeridian_geometry)
    
    # 2. Convert to Fiji's national metric grid (Fiji Red Beach, EPSG:2978)
    target_crs = "EPSG:2978"
    gdf = gdf.to_crs(target_crs)
    
    # 3. Calculate length in metres (now that we are in a metric CRS)
    gdf['length_m'] = gdf.geometry.length
    
    # 4. Ensure all are MultiLineString
    gdf['geometry'] = gdf['geometry'].apply(ensure_multilinestring)
    
    # 5. Save to GeoJSON
    gdf.to_file(output_file, driver='GeoJSON')
    print(f"Successfully processed and saved to {output_file}")

if __name__ == "__main__":
    process_transects('fiji_transects_wgs84.geojson', 'fiji_transects_fmg.geojson')