I am trying to figure out why I am getting inf for pyproj transformation from epsg:4326 to epsg:3857, I am using shapely to take in coordinates of a polygon so that I can rasterize. Afterwards, I rasterize it and get back points at which time I am trying to
transform using pyproj from epsg 4326 to epsg 3857. I have always_xy=True enabled to follow longitude, latitude formatting.

This is on very specific coordinates in Tuscon, AZ which so far as only been the only test case that resulted in inf value from pyproj.

I am wanting to expect results that are not inf

failed coordinates:

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

The coordiantes that work are:

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

Code:

`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}")`