Pyproj转换部分坐标返回inf值问题排查求助

我在使用Pyproj将EPSG:4326坐标转换为EPSG:3857时，部分亚利桑那州图森市的坐标会返回inf值，其他区域坐标转换正常。流程是先通过Shapely导入多边形坐标并栅格化，再对栅格化后的点做坐标转换，已经开启always_xy=True遵循经纬度格式，期望得到有效的转换结果。

失败坐标 #

[
{"lat": 32.1846842314981, "lng": -110.84070919725416},
{"lat": 32.184538948251245, "lng": -110.83227633211133},
{"lat": 32.17805544734469, "lng": -110.83210467073438},
{"lat": 32.177783021233786, "lng": -110.84083794328687}
]

正常坐标 #

[
{"lat": 40.6643340701724, "lng": -73.46681736057971},
{"lat": 40.65652089018049, "lng": -73.46630237644885},
{"lat": 40.65697668414809, "lng": -73.45540187901233},
{"lat": 40.66472470514906, "lng": -73.45617435520862}
]

原代码 #

import logging
from shapely.geometry import Polygon as ShapelyPolygon, Point
from pyproj import Transformer
from math import isfinite

# Setup logging
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

# Dummy settings with ZOOM_LEVEL
class settings:
    ZOOM_LEVEL = 18

# Define functions
def convert_to_shapely_polygon(coordinates):
    logger.info("Converting coordinates to Shapely polygon")
    try:
        polygon = ShapelyPolygon([(coord['lng'], coord['lat']) for coord in coordinates])
        logger.info(f"Converted polygon: {polygon}")
        return polygon
    except Exception as e:
        logger.error(f"Failed to convert coordinates to Shapely polygon. Error: {e}")
        raise Exception(f"Failed to convert coordinates to Shapely polygon: {str(e)}")

# Helper function to rasterize a polygon (Polygon)
def rasterize_polygon(polygon):
    try:
        minx, miny, maxx, maxy = polygon.bounds
        interval = .0027
        grid_points = []
        current_x = minx
        while current_x <= maxx:
            current_y = miny
            while current_y <= maxy:
                point = Point(current_x, current_y)
                if polygon.contains(point):
                    grid_points.append(point)
                current_y += interval
            current_x += interval

        logger.info(f"Rasterized polygon with {len(grid_points)} points")
        logger.info(f"Grid points: {grid_points}")
        return grid_points
    except Exception as e:
        logger.error(f"Failed to rasterize polygon. Error: {e}")
        raise Exception(f"Failed to rasterize polygon: {str(e)}")

def point_to_tile(point):
    try:
        transformer = Transformer.from_crs("epsg:4326", "epsg:3857", always_xy=True)
        logger.info(f"Transforming point: {point}")
        x_merc, y_merc = transformer.transform(point.y, point.x)

        logger.info(f"Transformed coordinates: x_merc={x_merc}, y_merc={y_merc}")

        # Add checks for infinity or NaN values
        if not (isfinite(x_merc) and isfinite(y_merc)):
            raise ValueError(f"Invalid transformed coordinates: x_merc={x_merc}, y_merc={y_merc}")

        n = 2.0 ** settings.ZOOM_LEVEL  # 2^zoom_level at 18
        tile_size_meters = 40075016.686 / n
        logger.info(f"tile_size_meters: {tile_size_meters}")

        x_tile = int((x_merc + 20037508.342) / tile_size_meters)
        y_tile = int((y_merc + 20037508.342) / tile_size_meters)
        logger.info(f"Tile coordinates: x_tile={x_tile}, y_tile={y_tile}")

        return x_tile, y_tile
    except OverflowError as e:
        logger.error(f"OverflowError: {e}")
        raise Exception("Error converting coordinates to tile: Overflow error")
    except Exception as e:
        logger.error(f"Unexpected error: {e}")
        raise Exception("Error converting coordinates to tile: Unexpected error")

# Helper function to convert tile coordinates to lat/lng (Polygon)
def tile_to_lat_lng(x_tile, y_tile):
    transformer = Transformer.from_crs("epsg:3857", "epsg:4326", always_xy=True)
    n = 2.0 ** settings.ZOOM_LEVEL  # 2^zoom_level at 18
    tile_size_meters = 40075016.686 / n
    x_merc = x_tile * tile_size_meters - 20037508.342
    y_merc = y_tile * tile_size_meters - 20037508.342
    lng, lat = transformer.transform(x_merc, y_merc)
    return lat, lng

# Main script logic
if __name__ == "__main__":
    coordinates_list = [
        [
            {"lat": 32.1846842314981, "lng": -110.84070919725416},
            {"lat": 32.184538948251245, "lng": -110.83227633211133},
            {"lat": 32.17805544734469, "lng": -110.83210467073438},
            {"lat": 32.177783021233786, "lng": -110.84083794328687}
        ],
        [
            {"lat": 40.6643340701724, "lng": -73.46681736057971},
            {"lat": 40.65652089018049, "lng": -73.46630237644885},
            {"lat": 40.65697668414809, "lng": -73.45540187901233},
            {"lat": 40.66472470514906, "lng": -73.45617435520862}
        ]
    ]

    try:
        for coordinates in coordinates_list:
            logger.info(f"Processing coordinates: {coordinates}")
            polygon = convert_to_shapely_polygon(coordinates)
            grid_points = rasterize_polygon(polygon)

            for point in grid_points:
                try:
                    x_tile, y_tile = point_to_tile(point)
                    print(f"Tile coordinates: x_tile={x_tile}, y_tile={y_tile}")
                    lat, lng = tile_to_lat_lng(x_tile, y_tile)
                    print(f"Lat/Lng coordinates: lat={lat}, lng={lng}")
                except Exception as e:
                    logger.error(f"Error processing point {point}: {e}")

    except Exception as e:
        logger.error(f"An error occurred: {e}")

问题原因及修复方案 #

核心问题 #

在point_to_tile函数中，调用transformer.transform时参数顺序错误：

• 开启always_xy=True后，Pyproj要求输入坐标顺序为**(经度, 纬度)**
• 原代码传入的是point.y, point.x，而Shapely的Point结构是Point(经度, 纬度)，即point.x是经度，point.y是纬度，参数顺序颠倒导致转换计算时出现数值溢出，返回inf。

修复后的关键代码 #

只需要修改point_to_tile函数中的转换参数行，同时修正Y瓦片坐标的计算逻辑（原逻辑会导致Y轴反转）：

def point_to_tile(point):
    try:
        transformer = Transformer.from_crs("epsg:4326", "epsg:3857", always_xy=True)
        logger.info(f"Transforming point: {point}")
        # 修正参数顺序：传入(point.x, point.y)即(经度, 纬度)
        x_merc, y_merc = transformer.transform(point.x, point.y)

        logger.info(f"Transformed coordinates: x_merc={x_merc}, y_merc={y_merc}")

        # Add checks for infinity or NaN values
        if not (isfinite(x_merc) and isfinite(y_merc)):
            raise ValueError(f"Invalid transformed coordinates: x_merc={x_merc}, y_merc={y_merc}")

        n = 2.0 ** settings.ZOOM_LEVEL  # 2^zoom_level at 18
        tile_size_meters = 40075016.686 / n
        logger.info(f"tile_size_meters: {tile_size_meters}")

        x_tile = int((x_merc + 20037508.342) / tile_size_meters)
        # 修正Y瓦片计算逻辑：Web墨卡托瓦片Y轴从北向南递减
        y_tile = int((20037508.342 - y_merc) / tile_size_meters)
        logger.info(f"Tile coordinates: x_tile={x_tile}, y_tile={y_tile}")

        return x_tile, y_tile
    except OverflowError as e:
        logger.error(f"OverflowError: {e}")
        raise Exception("Error converting coordinates to tile: Overflow error")
    except Exception as e:
        logger.error(f"Unexpected error: {e}")
        raise Exception("Error converting coordinates to tile: Unexpected error")

内容的提问来源于stack exchange，提问作者Jeremy Lin