0
# Coordinate Transformation
1

2
Functions for reprojecting and transforming map tiles and images between different coordinate reference systems. These capabilities enable integration with diverse geospatial workflows and coordinate systems.
3

4
## Capabilities
5

6
### Tile Warping
7

8
Reproject Web Mercator basemap tiles into any coordinate reference system with configurable resampling methods for optimal visual quality.
9

10
```python { .api }
11
def warp_tiles(img, extent, t_crs='EPSG:4326', resampling=Resampling.bilinear):
12
    """
13
    Reproject Web Mercator basemap into any CRS on-the-fly.
14

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

21
    Returns:
22
    tuple: (img: ndarray, extent: tuple)
23
        - img: Warped image as 3D array (height, width, bands)
24
        - extent: Bounding box (minX, maxX, minY, maxY) in target CRS
25
    """
26
```
27

28
**Usage Examples:**
29

30
```python
31
import contextily as ctx
32
from rasterio.enums import Resampling
33
import matplotlib.pyplot as plt
34

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

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

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

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

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

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

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

63
plt.show()
64
```
65

66
### Image Transform Warping
67

68
Reproject images with known transforms and coordinate systems, providing full control over the transformation process for advanced geospatial workflows.
69

70
```python { .api }
71
def warp_img_transform(img, transform, s_crs, t_crs, resampling=Resampling.bilinear):
72
    """
73
    Reproject image with given transform and source/target CRS.
74

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

82
    Returns:
83
    tuple: (w_img: ndarray, w_transform: affine.Affine)
84
        - w_img: Warped image as 3D array (bands, height, width)
85
        - w_transform: Transform of warped image
86
    """
87
```
88

89
**Usage Examples:**
90

91
```python
92
import contextily as ctx
93
import rasterio
94
from rasterio.transform import from_bounds
95
from rasterio.enums import Resampling
96
import numpy as np
97

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

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

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

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

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

123
# Save reprojected result
124
with rasterio.open('reprojected.tif', 'w',
125
                   driver='GTiff',
126
                   height=img_utm.shape[1],
127
                   width=img_utm.shape[2], 
128
                   count=img_utm.shape[0],
129
                   dtype=img_utm.dtype,
130
                   crs='EPSG:32618',
131
                   transform=transform_utm) as dst:
132
    dst.write(img_utm)
133
```
134

135
## Resampling Methods
136

137
### Available Resampling Algorithms
138

139
```python
140
from rasterio.enums import Resampling
141

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

151
### Choosing Resampling Methods
152

153
```python
154
import contextily as ctx
155
from rasterio.enums import Resampling
156

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

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

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

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

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

173
## Common Coordinate Systems
174

175
### Frequently Used CRS
176

177
```python
178
# Geographic coordinate systems (lat/lon)
179
'EPSG:4326'    # WGS84 - World wide standard
180
'EPSG:4269'    # NAD83 - North America  
181
'EPSG:4258'    # ETRS89 - Europe
182

183
# Projected coordinate systems  
184
'EPSG:3857'    # Web Mercator - Web mapping standard
185
'EPSG:32618'   # UTM Zone 18N - US East Coast
186
'EPSG:32633'   # UTM Zone 33N - Central Europe
187
'EPSG:3395'    # World Mercator
188
'EPSG:3031'    # Antarctic Polar Stereographic
189
'EPSG:3413'    # WGS 84 / NSIDC Sea Ice Polar Stereographic North
190
```
191

192
### CRS Format Options
193

194
```python
195
# String formats supported
196
ctx.warp_tiles(img, extent, t_crs='EPSG:4326')      # EPSG code
197
ctx.warp_tiles(img, extent, t_crs='WGS84')          # Common name
198
ctx.warp_tiles(img, extent, t_crs='+proj=longlat +datum=WGS84')  # PROJ string
199

200
# Using rasterio CRS objects
201
from rasterio.crs import CRS
202
target_crs = CRS.from_epsg(4326)
203
ctx.warp_tiles(img, extent, t_crs=target_crs)
204
```
205

206
## Integration Examples
207

208
### Working with GeoPandas
209

210
```python
211
import contextily as ctx
212
import geopandas as gpd
213
import matplotlib.pyplot as plt
214

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

219
# Get basemap in same CRS as data
220
if gdf.crs != 'EPSG:3857':
221
    # Get extent in data's CRS
222
    bounds = gdf.total_bounds  # [minx, miny, maxx, maxy]
223
    
224
    # Download tiles in Web Mercator
225
    img, extent_3857 = ctx.bounds2img(bounds[0], bounds[1], bounds[2], bounds[3], 
226
                                    zoom=12, ll=(gdf.crs == 'EPSG:4326'))
227
    
228
    # Transform tiles to match data CRS
229
    img_reproj, extent_reproj = ctx.warp_tiles(img, extent_3857, t_crs=gdf.crs)
230
    
231
    # Plot with matching coordinates
232
    ax = gdf.plot(figsize=(12, 8), alpha=0.7)
233
    ax.imshow(img_reproj, extent=extent_reproj, alpha=0.8)
234
    plt.show()
235
```
236

237
### Working with Rasterio
238

239
```python
240
import contextily as ctx
241
import rasterio
242
from rasterio.warp import transform_bounds
243

244
# Load existing raster data
245
with rasterio.open('analysis_area.tif') as src:
246
    data_crs = src.crs
247
    data_bounds = src.bounds
248
    
249
    # Transform bounds to WGS84 for tile download
250
    ll_bounds = transform_bounds(data_crs, 'EPSG:4326', *data_bounds)
251
    
252
    # Download basemap tiles
253
    img, extent_3857 = ctx.bounds2img(ll_bounds[0], ll_bounds[1], 
254
                                    ll_bounds[2], ll_bounds[3], 
255
                                    zoom=14, ll=True)
256
    
257
    # Transform tiles to match raster CRS
258
    img_reproj, extent_reproj = ctx.warp_tiles(img, extent_3857, t_crs=data_crs)
259
    
260
    # Now img_reproj and your raster data are in same CRS
261
    raster_data = src.read()
262
```
263

264
### Custom Projection Workflows
265

266
```python
267
import contextily as ctx
268
from pyproj import CRS, Transformer
269

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

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

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

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

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

288
# Web Mercator view
289
ax1.imshow(img, extent=extent)
290
ax1.set_title('Web Mercator')
291

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

296
plt.show()
297
```
298

299
## Performance Considerations
300

301
### Memory Management
302

303
```python
304
# Large image warping can be memory intensive
305
# Process in chunks for very large datasets
306
def warp_large_image(img, extent, t_crs, chunk_size=1000):
307
    """Warp large images in smaller chunks."""
308
    h, w, bands = img.shape
309
    warped_chunks = []
310
    
311
    for i in range(0, h, chunk_size):
312
        for j in range(0, w, chunk_size):
313
            # Extract chunk
314
            chunk = img[i:i+chunk_size, j:j+chunk_size, :]
315
            
316
            # Calculate chunk extent
317
            # ... extent calculation logic ...
318
            
319
            # Warp chunk
320
            chunk_warped, _ = ctx.warp_tiles(chunk, chunk_extent, t_crs=t_crs)
321
            warped_chunks.append(chunk_warped)
322
    
323
    # Reassemble chunks
324
    # ... reassembly logic ...
325
    return warped_image, warped_extent
326
```
327

328
### Quality vs Speed Tradeoffs
329

330
```python
331
from rasterio.enums import Resampling
332

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

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

341
# Balanced approach (good quality, reasonable speed)
342
img_balanced, extent_balanced = ctx.warp_tiles(img, extent,
343
                                             resampling=Resampling.cubic)
344
```
345

346
## Error Handling
347

348
### Common Transformation Issues
349

350
```python
351
import contextily as ctx
352
from rasterio.errors import CRSError
353

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

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

370
### Validation and Debugging
371

372
```python
373
import contextily as ctx
374
import rasterio
375
from rasterio.warp import transform_bounds
376

377
def validate_transformation(img, extent, source_crs, target_crs):
378
    """Validate transformation parameters before processing."""
379
    
380
    # Check if transformation is reasonable
381
    try:
382
        # Test transform with small extent
383
        test_bounds = transform_bounds(source_crs, target_crs, *extent)
384
        print(f"Source extent: {extent}")
385
        print(f"Target extent: {test_bounds}")
386
        
387
        # Check for extreme distortion
388
        source_area = (extent[1] - extent[0]) * (extent[3] - extent[2])
389
        target_area = (test_bounds[1] - test_bounds[0]) * (test_bounds[3] - test_bounds[2])
390
        area_ratio = target_area / source_area
391
        
392
        if area_ratio > 100 or area_ratio < 0.01:
393
            print(f"Warning: Large area distortion (ratio: {area_ratio:.2f})")
394
            
395
    except Exception as e:
396
        print(f"Transformation validation failed: {e}")
397
        return False
398
    
399
    return True
400

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