Coordinate Transformation

def warp_tiles(img, extent, t_crs='EPSG:4326', resampling=Resampling.bilinear):
    """
    Reproject Web Mercator basemap into any CRS on-the-fly.

    Parameters:
    - img: ndarray - Image as 3D array (height, width, bands) of RGB values
    - extent: tuple - Bounding box (minX, maxX, minY, maxY) in Web Mercator (EPSG:3857)  
    - t_crs: str or CRS - Target coordinate reference system (default: WGS84 lon/lat)
    - resampling: Resampling - Resampling method for warping (default: bilinear)

    Returns:
    tuple: (img: ndarray, extent: tuple)
        - img: Warped image as 3D array (height, width, bands)
        - extent: Bounding box (minX, maxX, minY, maxY) in target CRS
    """

import contextily as ctx
from rasterio.enums import Resampling
import matplotlib.pyplot as plt

# Download Web Mercator tiles
west, south, east, north = -74.1, 40.7, -73.9, 40.8  # NYC bounds (lon/lat)
img, extent_3857 = ctx.bounds2img(west, south, east, north, zoom=12, ll=True)

# Transform to WGS84 (lon/lat) coordinates
img_4326, extent_4326 = ctx.warp_tiles(img, extent_3857, t_crs='EPSG:4326')

# Transform to UTM Zone 18N (appropriate for NYC)
img_utm, extent_utm = ctx.warp_tiles(img, extent_3857, t_crs='EPSG:32618')

# Use different resampling method for better quality
img_cubic, extent_cubic = ctx.warp_tiles(img, extent_3857, 
                                       t_crs='EPSG:4326',
                                       resampling=Resampling.cubic)

# Plot comparison
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

ax1.imshow(img, extent=extent_3857)
ax1.set_title('Web Mercator (EPSG:3857)')
ax1.set_xlabel('Easting (m)')
ax1.set_ylabel('Northing (m)')

ax2.imshow(img_4326, extent=extent_4326)  
ax2.set_title('WGS84 (EPSG:4326)')
ax2.set_xlabel('Longitude')
ax2.set_ylabel('Latitude')

plt.show()

def warp_img_transform(img, transform, s_crs, t_crs, resampling=Resampling.bilinear):
    """
    Reproject image with given transform and source/target CRS.

    Parameters:
    - img: ndarray - Image as 3D array (bands, height, width) in rasterio format
    - transform: affine.Affine - Transform of input image from rasterio/affine
    - s_crs: str or CRS - Source coordinate reference system  
    - t_crs: str or CRS - Target coordinate reference system
    - resampling: Resampling - Resampling method for warping

    Returns:
    tuple: (w_img: ndarray, w_transform: affine.Affine)
        - w_img: Warped image as 3D array (bands, height, width)
        - w_transform: Transform of warped image
    """

import contextily as ctx
import rasterio
from rasterio.transform import from_bounds
from rasterio.enums import Resampling
import numpy as np

# Load existing georeferenced raster
with rasterio.open('input_raster.tif') as src:
    img = src.read()  # Shape: (bands, height, width)
    transform = src.transform
    source_crs = src.crs

# Reproject to different CRS
img_reproj, transform_reproj = ctx.warp_img_transform(
    img, transform, source_crs, 'EPSG:4326'
)

# Create transform from known bounds (Web Mercator example)
west_m, south_m, east_m, north_m = -8238000, 4969000, -8236000, 4971000
height, width = 1000, 1000
transform_wm = from_bounds(west_m, south_m, east_m, north_m, width, height)

# Create synthetic data for demonstration
synthetic_img = np.random.randint(0, 255, (3, height, width), dtype=np.uint8)

# Transform from Web Mercator to UTM
img_utm, transform_utm = ctx.warp_img_transform(
    synthetic_img, transform_wm, 'EPSG:3857', 'EPSG:32618',
    resampling=Resampling.nearest
)

# Save reprojected result
with rasterio.open('reprojected.tif', 'w',
                   driver='GTiff',
                   height=img_utm.shape[1],
                   width=img_utm.shape[2], 
                   count=img_utm.shape[0],
                   dtype=img_utm.dtype,
                   crs='EPSG:32618',
                   transform=transform_utm) as dst:
    dst.write(img_utm)

from rasterio.enums import Resampling

# Common resampling methods
Resampling.nearest     # Nearest neighbor - preserves pixel values, blocky appearance
Resampling.bilinear    # Bilinear interpolation - smooth, good for continuous data  
Resampling.cubic       # Cubic convolution - smoother, slower
Resampling.lanczos     # Lanczos windowed sinc - high quality, slower
Resampling.average     # Average of contributing pixels
Resampling.mode        # Most common value - good for categorical data

import contextily as ctx
from rasterio.enums import Resampling

# Download base image
img, extent = ctx.bounds2img(-74.1, 40.7, -73.9, 40.8, zoom=12, ll=True)

# Nearest neighbor - best for categorical/discrete data
img_nearest, _ = ctx.warp_tiles(img, extent, resampling=Resampling.nearest)

# Bilinear - good general purpose, default
img_bilinear, _ = ctx.warp_tiles(img, extent, resampling=Resampling.bilinear)

# Cubic - smoother for imagery, slower
img_cubic, _ = ctx.warp_tiles(img, extent, resampling=Resampling.cubic)

# Lanczos - highest quality for imagery
img_lanczos, _ = ctx.warp_tiles(img, extent, resampling=Resampling.lanczos)

# Geographic coordinate systems (lat/lon)
'EPSG:4326'    # WGS84 - World wide standard
'EPSG:4269'    # NAD83 - North America  
'EPSG:4258'    # ETRS89 - Europe

# Projected coordinate systems  
'EPSG:3857'    # Web Mercator - Web mapping standard
'EPSG:32618'   # UTM Zone 18N - US East Coast
'EPSG:32633'   # UTM Zone 33N - Central Europe
'EPSG:3395'    # World Mercator
'EPSG:3031'    # Antarctic Polar Stereographic
'EPSG:3413'    # WGS 84 / NSIDC Sea Ice Polar Stereographic North

# String formats supported
ctx.warp_tiles(img, extent, t_crs='EPSG:4326')      # EPSG code
ctx.warp_tiles(img, extent, t_crs='WGS84')          # Common name
ctx.warp_tiles(img, extent, t_crs='+proj=longlat +datum=WGS84')  # PROJ string

# Using rasterio CRS objects
from rasterio.crs import CRS
target_crs = CRS.from_epsg(4326)
ctx.warp_tiles(img, extent, t_crs=target_crs)

import contextily as ctx
import geopandas as gpd
import matplotlib.pyplot as plt

# Load vector data
gdf = gpd.read_file('data.shp')
print(f"Original CRS: {gdf.crs}")

# Get basemap in same CRS as data
if gdf.crs != 'EPSG:3857':
    # Get extent in data's CRS
    bounds = gdf.total_bounds  # [minx, miny, maxx, maxy]
    
    # Download tiles in Web Mercator
    img, extent_3857 = ctx.bounds2img(bounds[0], bounds[1], bounds[2], bounds[3], 
                                    zoom=12, ll=(gdf.crs == 'EPSG:4326'))
    
    # Transform tiles to match data CRS
    img_reproj, extent_reproj = ctx.warp_tiles(img, extent_3857, t_crs=gdf.crs)
    
    # Plot with matching coordinates
    ax = gdf.plot(figsize=(12, 8), alpha=0.7)
    ax.imshow(img_reproj, extent=extent_reproj, alpha=0.8)
    plt.show()

import contextily as ctx
import rasterio
from rasterio.warp import transform_bounds

# Load existing raster data
with rasterio.open('analysis_area.tif') as src:
    data_crs = src.crs
    data_bounds = src.bounds
    
    # Transform bounds to WGS84 for tile download
    ll_bounds = transform_bounds(data_crs, 'EPSG:4326', *data_bounds)
    
    # Download basemap tiles
    img, extent_3857 = ctx.bounds2img(ll_bounds[0], ll_bounds[1], 
                                    ll_bounds[2], ll_bounds[3], 
                                    zoom=14, ll=True)
    
    # Transform tiles to match raster CRS
    img_reproj, extent_reproj = ctx.warp_tiles(img, extent_3857, t_crs=data_crs)
    
    # Now img_reproj and your raster data are in same CRS
    raster_data = src.read()

import contextily as ctx
from pyproj import CRS, Transformer

# Define custom projection (example: Albers Equal Area for continental US)
custom_proj = """
+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 
+x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
"""

# Download tiles for continental US
img, extent = ctx.bounds2img(-125, 24, -66, 49, zoom=5, ll=True)

# Transform to custom projection
img_custom, extent_custom = ctx.warp_tiles(img, extent, t_crs=custom_proj)

# Use with projection-aware plotting
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))

# Web Mercator view
ax1.imshow(img, extent=extent)
ax1.set_title('Web Mercator')

# Custom projection view
ax2.imshow(img_custom, extent=extent_custom)
ax2.set_title('Albers Equal Area')

plt.show()

# Large image warping can be memory intensive
# Process in chunks for very large datasets
def warp_large_image(img, extent, t_crs, chunk_size=1000):
    """Warp large images in smaller chunks."""
    h, w, bands = img.shape
    warped_chunks = []
    
    for i in range(0, h, chunk_size):
        for j in range(0, w, chunk_size):
            # Extract chunk
            chunk = img[i:i+chunk_size, j:j+chunk_size, :]
            
            # Calculate chunk extent
            # ... extent calculation logic ...
            
            # Warp chunk
            chunk_warped, _ = ctx.warp_tiles(chunk, chunk_extent, t_crs=t_crs)
            warped_chunks.append(chunk_warped)
    
    # Reassemble chunks
    # ... reassembly logic ...
    return warped_image, warped_extent

from rasterio.enums import Resampling

# Fast transformation (lower quality)
img_fast, extent_fast = ctx.warp_tiles(img, extent, 
                                     resampling=Resampling.nearest)

# High quality transformation (slower)  
img_quality, extent_quality = ctx.warp_tiles(img, extent,
                                           resampling=Resampling.lanczos)

# Balanced approach (good quality, reasonable speed)
img_balanced, extent_balanced = ctx.warp_tiles(img, extent,
                                             resampling=Resampling.cubic)

import contextily as ctx
from rasterio.errors import CRSError

try:
    img_warped, extent_warped = ctx.warp_tiles(img, extent, t_crs='INVALID:1234')
except CRSError as e:
    print(f"Invalid CRS: {e}")
    # Use fallback CRS
    img_warped, extent_warped = ctx.warp_tiles(img, extent, t_crs='EPSG:4326')

# Handle extreme transformations
try:
    # Very large extent might cause issues
    img_warped, extent_warped = ctx.warp_tiles(img, extent, t_crs='EPSG:3413')
except Exception as e:
    print(f"Transformation failed: {e}")
    # Consider clipping to smaller extent or different resampling

import contextily as ctx
import rasterio
from rasterio.warp import transform_bounds

def validate_transformation(img, extent, source_crs, target_crs):
    """Validate transformation parameters before processing."""
    
    # Check if transformation is reasonable
    try:
        # Test transform with small extent
        test_bounds = transform_bounds(source_crs, target_crs, *extent)
        print(f"Source extent: {extent}")
        print(f"Target extent: {test_bounds}")
        
        # Check for extreme distortion
        source_area = (extent[1] - extent[0]) * (extent[3] - extent[2])
        target_area = (test_bounds[1] - test_bounds[0]) * (test_bounds[3] - test_bounds[2])
        area_ratio = target_area / source_area
        
        if area_ratio > 100 or area_ratio < 0.01:
            print(f"Warning: Large area distortion (ratio: {area_ratio:.2f})")
            
    except Exception as e:
        print(f"Transformation validation failed: {e}")
        return False
    
    return True

# Use validation before warping
if validate_transformation(img, extent, 'EPSG:3857', 'EPSG:4326'):
    img_warped, extent_warped = ctx.warp_tiles(img, extent, t_crs='EPSG:4326')