0
# Transformations
1

2
Functions for transforming raster images and vector data between different coordinate systems. Includes regridding, warping, mesh generation, and utilities for handling cyclic data.
3

4
## Capabilities
5

6
### Image Transformations
7

8
Transform raster images and arrays between different coordinate reference systems.
9

10
```python { .api }
11
def warp_array(array, target_proj, source_proj=None, target_res=(400, 200),
12
               source_extent=None, target_extent=None, mask_extrapolated=False):
13
    """
14
    Warp numpy array to target projection.
15
    
16
    Parameters:
17
    - array: numpy.ndarray, source data array
18
    - target_proj: cartopy.crs.CRS, target projection
19
    - source_proj: cartopy.crs.CRS, source projection (default PlateCarree)
20
    - target_res: tuple, target resolution (width, height)
21
    - source_extent: tuple, source extent (x0, x1, y0, y1)  
22
    - target_extent: tuple, target extent (x0, x1, y0, y1)
23
    - mask_extrapolated: bool, mask extrapolated values
24
    
25
    Returns:
26
    Tuple of (warped_array, target_extent)
27
    """
28

29
def warp_img(fname, target_proj, source_proj=None, target_res=(400, 200)):
30
    """
31
    Warp image file to target projection.
32
    
33
    Parameters:
34
    - fname: str, path to image file
35
    - target_proj: cartopy.crs.CRS, target projection
36
    - source_proj: cartopy.crs.CRS, source projection (default PlateCarree)
37
    - target_res: tuple, target resolution (width, height)
38
    
39
    Returns:
40
    Tuple of (warped_array, target_extent)
41
    """
42

43
def regrid(array, source_x_coords, source_y_coords, source_proj, target_proj,
44
           target_res, mask_extrapolated=False):
45
    """
46
    Regrid array between coordinate systems.
47
    
48
    Parameters:
49
    - array: numpy.ndarray, source data
50
    - source_x_coords: array-like, source x coordinates
51
    - source_y_coords: array-like, source y coordinates  
52
    - source_proj: cartopy.crs.CRS, source projection
53
    - target_proj: cartopy.crs.CRS, target projection
54
    - target_res: tuple, target resolution (width, height)
55
    - mask_extrapolated: bool, mask extrapolated values
56
    
57
    Returns:
58
    Tuple of (regridded_array, target_x_coords, target_y_coords)
59
    """
60

61
def mesh_projection(projection, nx, ny, x_extents=(None, None), y_extents=(None, None)):
62
    """
63
    Generate coordinate mesh for projection.
64
    
65
    Parameters:
66
    - projection: cartopy.crs.CRS, target projection
67
    - nx: int, number of x grid points
68
    - ny: int, number of y grid points
69
    - x_extents: tuple, x extent limits (min, max)
70
    - y_extents: tuple, y extent limits (min, max)
71
    
72
    Returns:
73
    Tuple of (x_mesh, y_mesh) coordinate arrays
74
    """
75
```
76

77
### Vector Field Transformations
78

79
Transform vector fields and scalar data to regular grids.
80

81
```python { .api }
82
def vector_scalar_to_grid(src_crs, target_proj, regrid_shape, x, y, u, v, *scalars, **kwargs):
83
    """
84
    Transform vector and scalar fields to regular grid.
85
    
86
    Parameters:
87
    - src_crs: cartopy.crs.CRS, source coordinate system
88
    - target_proj: cartopy.crs.CRS, target projection
89
    - regrid_shape: tuple, output grid shape (nx, ny)
90
    - x: array-like, source x coordinates
91
    - y: array-like, source y coordinates
92
    - u: array-like, x-component of vector field
93
    - v: array-like, y-component of vector field  
94
    - *scalars: additional scalar fields to transform
95
    - **kwargs: additional interpolation parameters
96
    
97
    Returns:
98
    Tuple of (target_x, target_y, target_u, target_v, *target_scalars)
99
    """
100
```
101

102
### Cyclic Data Utilities
103

104
Handle periodic/cyclic data boundaries common in global datasets.
105

106
```python { .api }
107
def add_cyclic_point(data, coord=None, axis=-1):
108
    """
109
    Add cyclic point to data array.
110
    
111
    Parameters:
112
    - data: array-like, input data array
113
    - coord: array-like, coordinate array (optional)
114
    - axis: int, axis along which to add cyclic point
115
    
116
    Returns:
117
    If coord provided: tuple of (extended_data, extended_coord)
118
    Otherwise: extended_data
119
    """
120

121
def add_cyclic(data, x=None, y=None, axis=-1, cyclic=360, precision=1e-4):
122
    """
123
    Add cyclic point with coordinate handling.
124
    
125
    Parameters:
126
    - data: array-like, input data
127
    - x: array-like, x coordinates (optional)
128
    - y: array-like, y coordinates (optional) 
129
    - axis: int, cyclic axis
130
    - cyclic: float, cyclic interval (default 360 for longitude)
131
    - precision: float, precision for cyclic detection
132
    
133
    Returns:
134
    Tuple of (extended_data, extended_x, extended_y) or subsets based on inputs
135
    """
136

137
def has_cyclic(x, axis=-1, cyclic=360, precision=1e-4):
138
    """
139
    Check if coordinates already have cyclic point.
140
    
141
    Parameters:
142
    - x: array-like, coordinate array
143
    - axis: int, axis to check
144
    - cyclic: float, expected cyclic interval
145
    - precision: float, precision for comparison
146
    
147
    Returns:
148
    bool, True if cyclic point exists
149
    """
150

151
def add_cyclic(data, x=None, y=None, axis=-1, cyclic=360, precision=1e-4):
152
    """
153
    Add cyclic point with coordinate handling (newer version of add_cyclic_point).
154
    
155
    Parameters:
156
    - data: array-like, input data
157
    - x: array-like, x coordinates (optional)
158
    - y: array-like, y coordinates (optional) 
159
    - axis: int, cyclic axis
160
    - cyclic: float, cyclic interval (default 360 for longitude)
161
    - precision: float, precision for cyclic detection
162
    
163
    Returns:
164
    Tuple of (extended_data, extended_x, extended_y) or subsets based on inputs
165
    """
166
```
167

168
### Geodesic Calculations
169

170
Geodesic computations on ellipsoidal Earth models.
171

172
```python { .api }
173
class Geodesic:
174
    """Geodesic calculations on ellipsoid."""
175
    def __init__(self, a, f):
176
        """
177
        Parameters:
178
        - a: float, semi-major axis
179
        - f: float, flattening
180
        """
181
    
182
    def inverse(self, lat1, lon1, lat2, lon2):
183
        """
184
        Compute inverse geodesic problem.
185
        
186
        Parameters:
187
        - lat1, lon1: float, first point coordinates
188
        - lat2, lon2: float, second point coordinates
189
        
190
        Returns:
191
        Dict with distance, azimuth angles
192
        """
193
    
194
    def direct(self, lat1, lon1, azi1, s12):
195
        """
196
        Compute direct geodesic problem.
197
        
198
        Parameters:
199
        - lat1, lon1: float, starting point
200
        - azi1: float, initial azimuth
201
        - s12: float, distance
202
        
203
        Returns:
204
        Dict with end point coordinates and final azimuth
205
        """
206
```
207

208
## Usage Examples
209

210
### Warping Images Between Projections
211

212
```python
213
import numpy as np
214
import matplotlib.pyplot as plt
215
import cartopy.crs as ccrs
216
from cartopy.img_transform import warp_array
217

218
# Create sample data on regular lat/lon grid
219
lons = np.linspace(-180, 180, 360)
220
lats = np.linspace(-90, 90, 180)
221
lon_grid, lat_grid = np.meshgrid(lons, lats)
222
data = np.sin(np.radians(lat_grid)) * np.cos(np.radians(lon_grid))
223

224
# Warp to different projections
225
source_proj = ccrs.PlateCarree()
226
target_proj = ccrs.Orthographic(central_longitude=-90, central_latitude=45)
227

228
warped_data, target_extent = warp_array(
229
    data, 
230
    target_proj, 
231
    source_proj=source_proj,
232
    target_res=(400, 400)
233
)
234

235
# Plot results
236
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
237

238
# Original data
239
ax1 = plt.subplot(1, 2, 1, projection=source_proj)
240
ax1.contourf(lon_grid, lat_grid, data, transform=source_proj)
241
ax1.coastlines()
242
ax1.set_global()
243
ax1.set_title('Original Data (PlateCarree)')
244

245
# Warped data  
246
ax2 = plt.subplot(1, 2, 2, projection=target_proj)
247
ax2.imshow(warped_data, extent=target_extent, transform=target_proj)
248
ax2.coastlines()
249
ax2.set_global()
250
ax2.set_title('Warped Data (Orthographic)')
251

252
plt.tight_layout()
253
plt.show()
254
```
255

256
### Handling Cyclic Data
257

258
```python
259
import numpy as np
260
import matplotlib.pyplot as plt
261
import cartopy.crs as ccrs
262
from cartopy.util import add_cyclic_point
263

264
# Create longitude data that doesn't include 360°
265
lons = np.arange(0, 360, 5)  # 0 to 355
266
lats = np.arange(-90, 91, 5)
267
lon_grid, lat_grid = np.meshgrid(lons, lats)
268

269
# Sample data
270
data = np.sin(np.radians(lat_grid)) * np.cos(np.radians(lon_grid))
271

272
# Add cyclic point to avoid gap at 180°/-180°
273
cyclic_data, cyclic_lons = add_cyclic_point(data, coord=lons)
274

275
# Plot comparison
276
fig, axes = plt.subplots(1, 2, figsize=(15, 6),
277
                        subplot_kw={'projection': ccrs.PlateCarree()})
278

279
# Without cyclic point
280
lon_mesh, lat_mesh = np.meshgrid(lons, lats)
281
axes[0].pcolormesh(lon_mesh, lat_mesh, data, transform=ccrs.PlateCarree())
282
axes[0].coastlines()
283
axes[0].set_global()
284
axes[0].set_title('Without Cyclic Point (Gap at Dateline)')
285

286
# With cyclic point
287
cyclic_lon_mesh, cyclic_lat_mesh = np.meshgrid(cyclic_lons, lats)
288
axes[1].pcolormesh(cyclic_lon_mesh, cyclic_lat_mesh, cyclic_data, 
289
                   transform=ccrs.PlateCarree())
290
axes[1].coastlines()
291
axes[1].set_global()
292
axes[1].set_title('With Cyclic Point (No Gap)')
293

294
plt.tight_layout()
295
plt.show()
296
```
297

298
### Vector Field Transformation
299

300
```python
301
import numpy as np
302
import matplotlib.pyplot as plt
303
import cartopy.crs as ccrs
304
from cartopy.vector_transform import vector_scalar_to_grid
305

306
# Create irregular vector field data
307
np.random.seed(42)
308
x_pts = np.random.uniform(-180, 180, 100)
309
y_pts = np.random.uniform(-90, 90, 100)
310
u_pts = 10 * np.sin(np.radians(y_pts))
311
v_pts = 5 * np.cos(np.radians(x_pts))
312

313
# Transform to regular grid
314
source_crs = ccrs.PlateCarree()
315
target_proj = ccrs.PlateCarree()
316

317
# Regular grid output
318
x_reg, y_reg, u_reg, v_reg = vector_scalar_to_grid(
319
    source_crs, target_proj, 
320
    regrid_shape=(20, 15),  # 20x15 regular grid
321
    x=x_pts, y=y_pts, u=u_pts, v=v_pts
322
)
323

324
# Plot results
325
fig = plt.figure(figsize=(15, 10))
326

327
# Original irregular data
328
ax1 = plt.subplot(2, 1, 1, projection=ccrs.PlateCarree())
329
ax1.quiver(x_pts, y_pts, u_pts, v_pts, transform=ccrs.PlateCarree(),
330
           alpha=0.7, color='blue', scale=200)
331
ax1.coastlines()
332
ax1.set_global()
333
ax1.set_title('Original Irregular Vector Field')
334

335
# Interpolated regular grid
336
ax2 = plt.subplot(2, 1, 2, projection=ccrs.PlateCarree())
337
ax2.quiver(x_reg, y_reg, u_reg, v_reg, transform=ccrs.PlateCarree(),
338
           color='red', scale=200)
339
ax2.coastlines()
340
ax2.set_global()
341
ax2.set_title('Interpolated Regular Grid')
342

343
plt.tight_layout()
344
plt.show()
345
```
346

347
### Creating Projection Meshes
348

349
```python
350
import numpy as np
351
import matplotlib.pyplot as plt
352
import cartopy.crs as ccrs
353
from cartopy.img_transform import mesh_projection
354

355
# Generate coordinate mesh for Lambert Conformal projection
356
proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=39)
357

358
# Create mesh over continental US region
359
x_mesh, y_mesh = mesh_projection(proj, nx=50, ny=40,
360
                                x_extents=(-2.5e6, 2.5e6),
361
                                y_extents=(-1.5e6, 1.5e6))
362

363
# Plot the mesh
364
fig = plt.figure(figsize=(12, 8))
365
ax = plt.axes(projection=proj)
366

367
# Plot mesh points
368
ax.scatter(x_mesh.flatten(), y_mesh.flatten(), s=1, 
369
           transform=proj, alpha=0.5)
370

371
# Add geographic features
372
ax.coastlines()
373
ax.add_feature(cfeature.BORDERS)
374
ax.set_extent([-125, -65, 25, 55], ccrs.PlateCarree())
375

376
plt.title('Lambert Conformal Projection Mesh')
377
plt.show()
378
```
379

380
### Geodesic Calculations
381

382
```python
383
from cartopy.geodesic import Geodesic
384
import matplotlib.pyplot as plt
385
import cartopy.crs as ccrs
386

387
# WGS84 ellipsoid parameters
388
geod = Geodesic(6378137.0, 1/298.257223563)
389

390
# Calculate great circle between two cities
391
# New York to London
392
ny_lat, ny_lon = 40.7128, -74.0060
393
london_lat, london_lon = 51.5074, -0.1278
394

395
result = geod.inverse(ny_lat, ny_lon, london_lat, london_lon)
396
distance_km = result['distance'] / 1000
397
azimuth = result['azi1']
398

399
print(f"Distance: {distance_km:.0f} km")
400
print(f"Initial bearing: {azimuth:.1f}°")
401

402
# Plot great circle path
403
fig = plt.figure(figsize=(12, 8))
404
ax = plt.axes(projection=ccrs.PlateCarree())
405

406
# Plot cities
407
ax.scatter([ny_lon, london_lon], [ny_lat, london_lat], 
408
           s=100, c=['red', 'blue'], transform=ccrs.PlateCarree())
409

410
# Add great circle (simplified - would need more points for accuracy)
411
ax.plot([ny_lon, london_lon], [ny_lat, london_lat], 
412
        'r--', transform=ccrs.Geodetic(), linewidth=2)
413

414
ax.coastlines()
415
ax.set_global()
416
ax.gridlines(draw_labels=True)
417
plt.title(f'Great Circle: New York to London ({distance_km:.0f} km)')
418
plt.show()
419
```
420

421
## Advanced Transformations
422

423
### Custom Interpolation
424

425
```python
426
from cartopy.vector_transform import vector_scalar_to_grid
427

428
# Custom interpolation with different methods
429
x_reg, y_reg, u_reg, v_reg, temp_reg = vector_scalar_to_grid(
430
    ccrs.PlateCarree(), ccrs.LambertConformal(),
431
    regrid_shape=(30, 25),
432
    x=lon_data, y=lat_data, u=u_wind, v=v_wind, temperature_data,
433
    method='linear',  # or 'nearest', 'cubic'
434
    fill_value=np.nan
435
)
436
```
437

438
### Masking Extrapolated Values
439

440
```python
441
# Warp with extrapolation masking
442
warped_data, extent = warp_array(
443
    original_data, target_projection,
444
    source_proj=ccrs.PlateCarree(),
445
    mask_extrapolated=True  # Mask values outside original domain
446
)
447

448
# Masked values will be np.nan
449
valid_mask = ~np.isnan(warped_data)
450
```