如何从GeoTIFF文件中提取海运航线的经纬度坐标？

Hey there! Sounds like you're already halfway there with the projection conversion to PlateCarree. Let's walk through how to pull out those actual shipping route coordinates from your 194M GeoTIFF file. We'll use Python's geospatial tools since you're already working with cartopy—super straightforward once you know the steps.

Step 1: Inspect the GeoTIFF Data #

First, we need to figure out which pixel values correspond to the shipping routes. Most shipping route GeoTIFFs use a specific value (like 1) to mark routes, with other values representing non-route areas. Let's read the file and check the unique pixel values:

import rasterio
import numpy as np

# Open your GeoTIFF file
with rasterio.open('raw_2008_shipping_mol_20150714093538.tif') as src:
    # Read the first band (adjust if your data uses a different band)
    route_data = src.read(1)
    # Confirm the projection is PlateCarree
    print(f"Dataset Projection: {src.crs}")
    # Grab the transform needed to convert pixels to coordinates
    pixel_to_geo = src.transform

# Check unique pixel values to find the one representing routes
unique_vals = np.unique(route_data)
print(f"Unique pixel values in the file: {unique_vals}")

Look at the output of unique_vals—pick the value that's associated with shipping routes (e.g., 1, 255, or another distinct number). Let's call this target_value for the next steps.

Step 2: Extract Route Pixels & Convert to Coordinates #

Now we'll find all pixels with your target value and convert their row/column positions to actual PlateCarree (latitude/longitude) coordinates:

Option 1: Extract Discrete Route Points #

If your GeoTIFF marks routes as individual points or dense pixel clusters, this method works:

from rasterio.transform import xy

# Find all row/column indices where the pixel matches the route value
row_idx, col_idx = np.where(route_data == target_value)

# Convert pixel indices to lat/lon coordinates
longitudes = []
latitudes = []
for row, col in zip(row_idx, col_idx):
    lon, lat = xy(pixel_to_geo, row, col)
    longitudes.append(lon)
    latitudes.append(lat)

# Optional: Organize into a DataFrame for easier handling
import pandas as pd
route_coords = pd.DataFrame({"Longitude": longitudes, "Latitude": latitudes})
print(route_coords.head())

Option 2: Extract Continuous Route Lines/Polygons #

If your routes are stored as continuous linear regions (not just scattered points), using vectorization will give you cleaner, connected route geometries:

from rasterio.features import shapes
import geopandas as gpd

# Extract vector shapes from the raster data
route_shapes = (
    {"geometry": shape, "properties": {"value": val}}
    for shape, val in shapes(route_data, transform=pixel_to_geo)
    if val == target_value
)

# Convert to a GeoDataFrame (geospatial-enabled DataFrame)
route_gdf = gpd.GeoDataFrame.from_features(list(route_shapes))
# Set the coordinate reference system to PlateCarree
route_gdf.crs = src.crs

# Extract lat/lon from the geometry column if needed
route_gdf["Longitude"] = route_gdf.geometry.x
route_gdf["Latitude"] = route_gdf.geometry.y

print(route_gdf.head())

Quick Notes to Troubleshoot #

• Projection Check: If your initial projection conversion was done separately, make sure the GeoTIFF you're reading is already in cartopy.PlateCarree() (EPSG:4326). If not, you can reproject it using rasterio.warp.reproject first.
• Band Selection: Some GeoTIFFs use multiple bands—double-check that you're reading the correct band with the route data.
• Large File Handling: Since your file is 194M, using numpy's vectorized operations (like np.where) will be faster than looping through every pixel manually.

Let me know if you hit snags with identifying the target pixel value or dealing with projection mismatches—I can help refine this further!

内容的提问来源于stack exchange，提问作者Li Ziming