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("CRS:", gdf.crs)
    print("Geometry types:", gdf.geometry.type.unique())
    print("First 5 rows:")
    print(gdf.head())
    print("\nSummary statistics for geometry:")
    print(gdf.geometry.describe())

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

import geopandas as gpd

def inspect_geojson(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("CRS:", gdf.crs)
    print("Geometry types:", gdf.geometry.type.unique())
    print("First 5 rows:")
    print(gdf.head())
    print("\nSummary statistics for geometry:")
    print(gdf.geometry.describe())

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

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

# Fiji National Grid: EPSG:2977 (Fiji Grid)
# Wait, let me check if EPSG:2977 is correct or if there's a better one.
# Actually, let's use a common projection for Fiji if 2977 is not available or search.
# Commonly, Fiji uses UTM zones. For Fiji, UTM zones 58S and 59S cover much of it.
# But the prompt specifies "Fiji's national metric grid". 
# Let's assume EPSG:2977 is the Fiji Grid.

def fix_and_convert(input_file, output_file):
    gdf = gpd.read_file(input_file)
    
    # The problem: lines spanning the globe.
    # This often happens with longitudes crossing the antimeridian (180/-180).
    # If a LineString goes from 179 to -179, it's drawn as a long arc across the globe.
    # We need to split these at the antimeridian.
    
    new_geometries = []
    
    for geom in gdf.geometry:
        if geom.is_empty:
            new_geometries.append(geom)
            continue
            
        coords = list(geom.coords)
        split_segments = []
        current_segment = [coords[0]]
        
        for i in range(1, len(coords)):
            prev_lon = coords[i-1][0]
            curr_lon = coords[i][0]
            
            # Check for antimeridian crossing. 
            # A jump of more than 180 degrees in longitude indicates a wrap-around.
            if abs(curr_lon - prev_lon) > 180:
                # Normalize longitudes to be continuous
                # If prev is 179 and curr is -179, curr should be treated as 181
                # or prev as -181. 
                # However, a simpler way to represent this for MultiLineString:
                # split the segment at the antimeridian.
                
                # Let's try a more robust approach: 
                # If we detect a jump, we wrap the coordinate.
                # But the requirement says "fix the geometries so they render correctly".
                # The problem is the LINESTRING goes the "long way" around the world.
                # We want the "short way".
                
                # Let's implement the split.
                # We'll create a segment from current_segment to the point before the jump,
                # and start a new segment.
                
                # But we need to handle the longitude wrap to make the segments "local".
                # If the jump is from 179 to -179, the segment should be 179 to 180 (effectively),
                # then -180 to -179.
                
                # Actually, for MultiLineString, we can just split it.
                # To avoid the "long arc", we need to ensure the coordinates are adjusted 
                # so they don't jump. But if we adjust them, they are no longer WGS84.
                # Wait, if we split it, we'll have two segments:
                # One ending at 180 and one starting at -180.
                # But WGS84 is bounded -180 to 180.
                
                # Actually, standard way to handle this in GIS is to split the line at the antimeridian.
                # Segment 1: ... -> 180
                # Segment 2: -180 -> ...
                
                # Let's find the exact point of crossing.
                # For simplicity in this task, if we detect a jump > 180, 
                # we split the LineString at the index of the jump.
                
                # However, we must also ensure the jump doesn't result in an "empty" 
                # segment if we just split.
                
                # Let's refine:
                # If |lon2 - lon1| > 180:
                #   it's an antimeridian crossing.
                #   We split the LineString into two (or more) parts.
                #   But wait, if the user wants them as MultiLineString, 
                #   we can just split the LineString at the jump.
                #   The resulting segments will be in the -180/180 range.
                
                # Example: (179, -17) to (-179, -17)
                # Segment 1: (179, -17) to (180, -17) ?? No, that's not quite right.
                # The actual geometry provided has a jump.
                # We can split it at the index where the jump occurs.
                # But the coordinates themselves in the file might be -179.
                # If we just split, we get (179, -17) and (-179, -17).
                # This is correct and doesn't form a long arc when rendered in a 
                # system that understands the antimeridian.
                # BUT the prompt says "some of the lines... render as arcs spanning the entire globe".
                # This usually happens because the renderer follows the line from 179 to -179.
                # By splitting it into two LineStrings, we prevent that arc.
                
                split_segments.append(LineString(current_segment))
                current_segment = [coords[i]]
            else:
                current_segment.append(coords[i])
        
        if current_segment:
            split_segments.append(LineString(current_segment))
            
        if len(split_segments) > 1:
            new_geometries.append(MultiLineString(split_segments))
        elif len(split_segments) == 1:
            new_geometries.append(split_segments[0])
        else:
            new_geometries.append(geom) # Should not happen
            
    gdf.geometry = new_geometries
    
    # Now convert to Fiji's national metric grid.
    # I'll use EPSG:2977 (Fiji Grid) if available, otherwise I'll check for a standard one.
    # Let's try to project. 
    # If EPSG:2977 fails, I'll use a UTM zone appropriate for Fiji.
    # Fiji is around 177-180E, 17-21S.
    # UTM Zone 58S or 59S.
    
    try:
        gdf_projected = gdf.to_crs(epsg=2977)
        print("Successfully projected to EPSG:2977")
    except Exception as e:
        print(f"EPSG:2977 failed: {e}. Trying UTM Zone 58S (EPSG:32758)...")
        gdf_projected = gdf.to_crs(epsg=32758)
        print("Successfully projected to EPSG:32758")

    # Calculate length in meters
    gdf_projected['length_m'] = gdf_projected.geometry.length
    
    # Ensure all geometries are MultiLineString
    def ensure_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    gdf_projected['geometry'] = gdf_projected.geometry.apply(ensure_multilinestring)
    
    # Final cleanup: ensure we only keep required columns and order them? 
    # The prompt says "Carry the survey attributes (vessel and survey_date) through unchanged, 
    # use transect_id as the key field". It doesn't say EXCLUDE others, 
    # but it's good practice to keep them.
    
    # Write to GeoJSON
    gdf_projected.to_file(output_file, driver='GeoJSON')
    print(f"Saved to {output_file}")

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

import geopandas as gpd

def verify_output(file_path):
    gdf = gpd.read_file(file_path)
    print("Columns:", gdf.columns.tolist())
    print("CRS:", gdf.crs)
    print("Geometry types:", gdf.geometry.type.unique())
    print("First 5 rows:")
    print(gdf.head())
    print("\nCheck for length_m column:")
    print("length_m" in gdf.columns)
    print("All geometries are MultiLineString?")
    print(all(geom.geom_type == 'MultiLineString' for geom in gdf.geometry))

if __name__ == "__main__":
    verify_output("fiji_transects_fmg.geojson")

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

def fix_and_convert(input_file, output_file):
    gdf = gpd.read_file(input_file)
    
    new_geometries = []
    
    for geom in gdf.geometry:
        if geom.is_empty:
            new_geometries.append(geom)
            continue
            
        coords = list(geom.coords)
        split_segments = []
        current_segment = [coords[0]]
        
        for i in range(1, len(coords)):
            prev_lon = coords[i-1][0]
            curr_lon = coords[i][0]
            
            # Antimeridian crossing detection: 
            # If the longitude jump is greater than 180, it's a wrap-around.
            if abs(curr_lon - prev_lon) > 180:
                if len(current_segment) > 1:
                    split_segments.append(LineString(current_segment))
                current_segment = [coords[i]]
            else:
                current_segment.append(coords[i])
        
        if len(current_segment) > 1:
            split_segments.append(LineString(current_segment))
        elif len(current_segment) == 1:
            # A single point shouldn't happen for a LineString, but for safety:
            pass

        if not split_segments:
            # If for some reason no segments were created (shouldn't happen for non-empty)
            new_geometries.append(geom)
        elif len(split_segments) == 1:
            # We still want to convert it to MultiLineString as per requirement
            new_geometries.append(MultiLineString([split_segments[0]]))
        else:
            new_geometries.append(MultiLineString(split_segments))
            
    gdf.geometry = new_geometries
    
    # Project to Fiji's national metric grid (EPSG:2977)
    gdf_projected = gdf.to_crs(epsg=2977)

    # Calculate length in meters
    gdf_projected['length_m'] = gdf_projected.geometry.length
    
    # Final check to ensure every feature is MultiLineString
    def ensure_multilinestring(geom):
        if geom.geom_type == 'LineString':
            return MultiLineString([geom])
        return geom

    gdf_projected['geometry'] = gdf_projected.geometry.apply(ensure_multilinestring)
    
    # Save the result
    gdf_projected.to_file(output_file, driver='GeoJSON')

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