Visualization of map_raster test

This notebook visualizes the test for map_raster using: - Real S1A SAR coordinates (Gulf of California) - Synthetic ECMWF wind data

Source: tests/test_map_raster_4cases.py (comprehensive 4-case test suite)

import sys
from pathlib import Path

# Add parent directories to path to import from tests and mapraster
sys.path.insert(0, str(Path.cwd().parent.parent))

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from mapraster.main import map_raster
from tests.tools_test import fake_ecmwf_0100_1h, build_footprint

%matplotlib inline

1. Create SAR dataset with real S1A coordinates

def create_real_sar_dataset():
    """
    Create SAR dataset with real S1A coordinates from Gulf of California.
    Based on s1a-iw-owi-dv-20210909t130650-20210909t130715-039605-04AE83.nc

    Returns synthetic dataset with bilinear interpolation between corner points.
    """
    # Real corner coordinates from S1A SAR
    corners = {
        'lon': np.array([-106.73246, -109.26672, -109.43854, -106.90172]),
        'lat': np.array([21.72079, 22.14191, 20.64836, 20.22559])
    }

    # SAR grid dimensions
    nlines, nsamples = 167, 256

    # Create meshgrid indices
    lines = np.arange(nlines)
    samples = np.arange(nsamples)

    # Bilinear interpolation
    s_norm = samples / (nsamples - 1)
    l_norm = lines / (nlines - 1)

    # Top edge (line=0): interpolate between corners[0] and corners[1]
    lon_top = corners['lon'][0] + s_norm * (corners['lon'][1] - corners['lon'][0])
    lat_top = corners['lat'][0] + s_norm * (corners['lat'][1] - corners['lat'][0])

    # Bottom edge (line=-1): interpolate between corners[3] and corners[2]
    lon_bottom = corners['lon'][3] + s_norm * (corners['lon'][2] - corners['lon'][3])
    lat_bottom = corners['lat'][3] + s_norm * (corners['lat'][2] - corners['lat'][3])

    # Full grid: interpolate between top and bottom edges
    lon_grid = lon_top[None, :] + l_norm[:, None] * (lon_bottom[None, :] - lon_top[None, :])
    lat_grid = lat_top[None, :] + l_norm[:, None] * (lat_bottom[None, :] - lat_top[None, :])

    # Create xarray Dataset
    ds = xr.Dataset(
        {
            'longitude': (['line', 'sample'], lon_grid),
            'latitude': (['line', 'sample'], lat_grid),
        },
        coords={
            'line': lines,
            'sample': samples,
        }
    )

    return ds

sar_dataset = create_real_sar_dataset()
print(f"SAR shape: {sar_dataset['longitude'].shape}")
print(f"SAR lon range: [{sar_dataset['longitude'].min().values:.2f}, {sar_dataset['longitude'].max().values:.2f}]")
print(f"SAR lat range: [{sar_dataset['latitude'].min().values:.2f}, {sar_dataset['latitude'].max().values:.2f}]")

SAR shape: (167, 256)
SAR lon range: [-109.44, -106.73]
SAR lat range: [20.23, 22.14]

sar_dataset

<xarray.Dataset> Size: 687kB
Dimensions:    (line: 167, sample: 256)
Coordinates:
  * line       (line) int64 1kB 0 1 2 3 4 5 6 7 ... 160 161 162 163 164 165 166
  * sample     (sample) int64 2kB 0 1 2 3 4 5 6 ... 249 250 251 252 253 254 255
Data variables:
    longitude  (line, sample) float64 342kB -106.7 -106.7 ... -109.4 -109.4
    latitude   (line, sample) float64 342kB 21.72 21.72 21.72 ... 20.65 20.65

xarray.Dataset

• Dimensions:
  • line: 167
  • sample: 256
• Coordinates: (2)
  • line
    (line)
    int64
    0 1 2 3 4 5 ... 162 163 164 165 166

    array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,
            14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
            28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,  40,  41,
            42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  55,
            56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  66,  67,  68,  69,
            70,  71,  72,  73,  74,  75,  76,  77,  78,  79,  80,  81,  82,  83,
            84,  85,  86,  87,  88,  89,  90,  91,  92,  93,  94,  95,  96,  97,
            98,  99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111,
           112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
           126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139,
           140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153,
           154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166])

  • sample
    (sample)
    int64
    0 1 2 3 4 5 ... 251 252 253 254 255

    array([  0,   1,   2, ..., 253, 254, 255], shape=(256,))

• Data variables: (2)
  • longitude
    (line, sample)
    float64
    -106.7 -106.7 ... -109.4 -109.4

    array([[-106.73246   , -106.74239827, -106.75233655, ..., -109.24684345,
            -109.25678173, -109.26672   ],
           [-106.73347964, -106.74341797, -106.75335631, ..., -109.24787839,
            -109.25781673, -109.26775506],
           [-106.73449928, -106.74443767, -106.75437607, ..., -109.24891333,
            -109.25885173, -109.26879012],
           ...,
           [-106.89968072, -106.90962892, -106.91957711, ..., -109.41657349,
            -109.42652169, -109.43646988],
           [-106.90070036, -106.91064861, -106.92059687, ..., -109.41760843,
            -109.42755669, -109.43750494],
           [-106.90172   , -106.91166831, -106.92161663, ..., -109.41864337,
            -109.42859169, -109.43854   ]], shape=(167, 256))

  • latitude
    (line, sample)
    float64
    21.72 21.72 21.72 ... 20.65 20.65

    array([[21.72079   , 21.72244145, 21.7240929 , ..., 22.1386071 ,
            22.14025855, 22.14191   ],
           [21.71178277, 21.71343426, 21.71508575, ..., 22.12960973,
            22.13126122, 22.13291271],
           [21.70277554, 21.70442707, 21.7060786 , ..., 22.12061236,
            22.12226389, 22.12391542],
           ...,
           [20.24360446, 20.2452623 , 20.24692015, ..., 20.66303889,
            20.66469673, 20.66635458],
           [20.23459723, 20.23625511, 20.23791299, ..., 20.65404152,
            20.65569941, 20.65735729],
           [20.22559   , 20.22724792, 20.22890584, ..., 20.64504416,
            20.64670208, 20.64836   ]], shape=(167, 256))

2. Create synthetic ECMWF wind data

ecmwf_dataset = fake_ecmwf_0100_1h(to180=True, with_nan=False)
print(f"ECMWF shape: {ecmwf_dataset['U10'].shape}")
print(f"ECMWF x range: [{ecmwf_dataset['x'].min().values:.2f}, {ecmwf_dataset['x'].max().values:.2f}]")
print(f"ECMWF y range: [{ecmwf_dataset['y'].min().values:.2f}, {ecmwf_dataset['y'].max().values:.2f}]")

ECMWF shape: (1810, 3600)
ECMWF x range: [-180.00, 179.90]
ECMWF y range: [-90.00, 90.00]

ecmwf_dataset

<xarray.Dataset> Size: 104MB
Dimensions:      (y: 1810, x: 3600)
Coordinates:
  * y            (y) float64 14kB -90.0 -89.9 -89.8 -89.7 ... 89.8 89.9 90.0
  * x            (x) float64 29kB -180.0 -179.9 -179.8 ... 179.7 179.8 179.9
    spatial_ref  int64 8B 0
Data variables:
    U10          (y, x) float64 52MB 5.0 5.0 5.0 5.0 5.0 ... 5.0 5.0 5.0 5.0 5.0
    V10          (y, x) float64 52MB 2.449e-16 -0.003491 ... 0.006981 0.003491
Attributes:
    time:     2023-03-03 07:00:00

xarray.Dataset

• Dimensions:
  • y: 1810
  • x: 3600
• Coordinates: (3)
  • y
    (y)
    float64
    -90.0 -89.9 -89.8 ... 89.9 90.0

    array([-90.      , -89.900498, -89.800995, ...,  89.800995,  89.900498,
            90.      ], shape=(1810,))

  • x
    (x)
    float64
    -180.0 -179.9 ... 179.8 179.9

    array([-180. , -179.9, -179.8, ...,  179.7,  179.8,  179.9], shape=(3600,))

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Data variables: (2)
  • U10
    (y, x)
    float64
    5.0 5.0 5.0 5.0 ... 5.0 5.0 5.0 5.0

    array([[5.        , 5.        , 5.        , ..., 5.        , 5.        ,
            5.        ],
           [5.00347329, 5.00347329, 5.00347329, ..., 5.00347329, 5.00347329,
            5.00347329],
           [5.00694657, 5.00694657, 5.00694657, ..., 5.00694657, 5.00694657,
            5.00694657],
           ...,
           [5.00694657, 5.00694657, 5.00694657, ..., 5.00694657, 5.00694657,
            5.00694657],
           [5.00347329, 5.00347329, 5.00347329, ..., 5.00347329, 5.00347329,
            5.00347329],
           [5.        , 5.        , 5.        , ..., 5.        , 5.        ,
            5.        ]], shape=(1810, 3600))

  • V10
    (y, x)
    float64
    2.449e-16 -0.003491 ... 0.003491

    array([[ 2.44929360e-16, -3.49065673e-03, -6.98130283e-03, ...,
             1.04719277e-02,  6.98130283e-03,  3.49065673e-03],
           [ 2.44929360e-16, -3.49065673e-03, -6.98130283e-03, ...,
             1.04719277e-02,  6.98130283e-03,  3.49065673e-03],
           [ 2.44929360e-16, -3.49065673e-03, -6.98130283e-03, ...,
             1.04719277e-02,  6.98130283e-03,  3.49065673e-03],
           ...,
           [ 2.44929360e-16, -3.49065673e-03, -6.98130283e-03, ...,
             1.04719277e-02,  6.98130283e-03,  3.49065673e-03],
           [ 2.44929360e-16, -3.49065673e-03, -6.98130283e-03, ...,
             1.04719277e-02,  6.98130283e-03,  3.49065673e-03],
           [ 2.44929360e-16, -3.49065673e-03, -6.98130283e-03, ...,
             1.04719277e-02,  6.98130283e-03,  3.49065673e-03]],
          shape=(1810, 3600))

• Attributes: (1)

  time :

  2023-03-03 07:00:00

3. Run map_raster

footprint = build_footprint(sar_dataset)

result = map_raster(
    raster_ds=ecmwf_dataset,
    originalDataset=sar_dataset,
    footprint=footprint,
    cross_antimeridian=False,
)

print(f"Result shape: {result['U10'].shape}")
print(f"Result U10 range: [{result['U10'].min().values:.2f}, {result['U10'].max().values:.2f}]")
print(f"Result V10 range: [{result['V10'].min().values:.2f}, {result['V10'].max().values:.2f}]")

Result shape: (167, 256)
Result U10 range: [6.85, 6.88]
Result V10 range: [-1.92, -1.89]

4. Visualizations

4.1 SAR grid geometry

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Longitude
im1 = ax1.pcolormesh(
    sar_dataset['sample'],
    sar_dataset['line'],
    sar_dataset['longitude'],
    shading='auto',
    cmap='viridis'
)
ax1.set_xlabel('Sample')
ax1.set_ylabel('Line')
ax1.set_title('SAR Longitude Grid')
plt.colorbar(im1, ax=ax1, label='Longitude (°)')

# Latitude
im2 = ax2.pcolormesh(
    sar_dataset['sample'],
    sar_dataset['line'],
    sar_dataset['latitude'],
    shading='auto',
    cmap='plasma'
)
ax2.set_xlabel('Sample')
ax2.set_ylabel('Line')
ax2.set_title('SAR Latitude Grid')
plt.colorbar(im2, ax=ax2, label='Latitude (°)')

plt.tight_layout()
plt.show()

4.2 ECMWF wind fields (global)

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(14, 12))

# Wind speed
wind_speed_ecmwf = np.sqrt(ecmwf_dataset['U10']**2 + ecmwf_dataset['V10']**2)
im0 = ax1.pcolormesh(
    ecmwf_dataset['x'],
    ecmwf_dataset['y'],
    wind_speed_ecmwf,
    shading='auto',
    cmap='viridis',
    vmin=0,
    vmax=15
)
ax1.set_xlabel('Longitude (°)')
ax1.set_ylabel('Latitude (°)')
ax1.set_title('ECMWF Wind Speed (global)')
plt.colorbar(im0, ax=ax1, label='Wind Speed (m/s)')

# U10
im1 = ax2.pcolormesh(
    ecmwf_dataset['x'],
    ecmwf_dataset['y'],
    ecmwf_dataset['U10'],
    shading='auto',
    cmap='RdBu_r',
    vmin=-10,
    vmax=10
)
ax2.set_xlabel('Longitude (°)')
ax2.set_ylabel('Latitude (°)')
ax2.set_title('ECMWF U10 Wind Component (global)')
plt.colorbar(im1, ax=ax2, label='U10 (m/s)')

# V10
im2 = ax3.pcolormesh(
    ecmwf_dataset['x'],
    ecmwf_dataset['y'],
    ecmwf_dataset['V10'],
    shading='auto',
    cmap='RdBu_r',
    vmin=-10,
    vmax=10
)
ax3.set_xlabel('Longitude (°)')
ax3.set_ylabel('Latitude (°)')
ax3.set_title('ECMWF V10 Wind Component (global)')
plt.colorbar(im2, ax=ax3, label='V10 (m/s)')

plt.tight_layout()
plt.show()

4.3 ECMWF wind fields (SAR region zoom)

# Extract SAR region from ECMWF
lon_min, lon_max = sar_dataset['longitude'].min().values, sar_dataset['longitude'].max().values
lat_min, lat_max = sar_dataset['latitude'].min().values, sar_dataset['latitude'].max().values

ecmwf_subset = ecmwf_dataset.sel(
    x=slice(lon_min-1, lon_max+1),
    y=slice(lat_min-1, lat_max+1)
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# U10 zoom
im1 = ax1.pcolormesh(
    ecmwf_subset['x'],
    ecmwf_subset['y'],
    ecmwf_subset['U10'],
    shading='auto',
    cmap='RdBu_r'
)
ax1.set_xlabel('Longitude (°)')
ax1.set_ylabel('Latitude (°)')
ax1.set_title('ECMWF U10 (SAR region)')
plt.colorbar(im1, ax=ax1, label='U10 (m/s)')

# Add SAR footprint
lon_corners = sar_dataset['longitude'].values[[0, 0, -1, -1, 0], [0, -1, -1, 0, 0]]
lat_corners = sar_dataset['latitude'].values[[0, 0, -1, -1, 0], [0, -1, -1, 0, 0]]
ax1.plot(lon_corners, lat_corners, 'k-', linewidth=2, label='SAR footprint')
ax1.legend()

# V10 zoom
im2 = ax2.pcolormesh(
    ecmwf_subset['x'],
    ecmwf_subset['y'],
    ecmwf_subset['V10'],
    shading='auto',
    cmap='RdBu_r'
)
ax2.set_xlabel('Longitude (°)')
ax2.set_ylabel('Latitude (°)')
ax2.set_title('ECMWF V10 (SAR region)')
plt.colorbar(im2, ax=ax2, label='V10 (m/s)')

# Add SAR footprint
ax2.plot(lon_corners, lat_corners, 'k-', linewidth=2, label='SAR footprint')
ax2.legend()

plt.tight_layout()
plt.show()

4.4 Interpolated wind on SAR grid (result of map_raster)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# U10 interpolated
im1 = ax1.pcolormesh(
    sar_dataset['sample'],
    sar_dataset['line'],
    result['U10'],
    shading='auto',
    cmap='RdBu_r'
)
ax1.set_xlabel('Sample')
ax1.set_ylabel('Line')
ax1.set_title('Interpolated U10 on SAR grid')
plt.colorbar(im1, ax=ax1, label='U10 (m/s)')

# V10 interpolated
im2 = ax2.pcolormesh(
    sar_dataset['sample'],
    sar_dataset['line'],
    result['V10'],
    shading='auto',
    cmap='RdBu_r'
)
ax2.set_xlabel('Sample')
ax2.set_ylabel('Line')
ax2.set_title('Interpolated V10 on SAR grid')
plt.colorbar(im2, ax=ax2, label='V10 (m/s)')

plt.tight_layout()
plt.show()

4.5 Before/After comparison in lon/lat coordinates

fig, axes = plt.subplots(2, 2, figsize=(16, 12))

# Extract SAR coordinates
sar_lon = sar_dataset['longitude'].values
sar_lat = sar_dataset['latitude'].values

# U10: Before (ECMWF subset)
ax1 = axes[0, 0]
im1 = ax1.pcolormesh(
    ecmwf_subset['x'],
    ecmwf_subset['y'],
    ecmwf_subset['U10'],
    shading='auto',
    cmap='RdBu_r'
)
ax1.plot(lon_corners, lat_corners, 'k-', linewidth=2, label='SAR footprint')
ax1.set_xlabel('Longitude (°)')
ax1.set_ylabel('Latitude (°)')
ax1.set_title('ECMWF U10 (Before interpolation)')
ax1.legend()
plt.colorbar(im1, ax=ax1, label='U10 (m/s)')

# U10: After (interpolated on SAR grid)
ax2 = axes[0, 1]
im2 = ax2.scatter(
    sar_lon.flatten(),
    sar_lat.flatten(),
    c=result['U10'].values.flatten(),
    s=5,
    cmap='RdBu_r',
    vmin=im1.get_clim()[0],
    vmax=im1.get_clim()[1]
)
ax2.set_xlabel('Longitude (°)')
ax2.set_ylabel('Latitude (°)')
ax2.set_title('Interpolated U10 (After map_raster)')
plt.colorbar(im2, ax=ax2, label='U10 (m/s)')

# V10: Before (ECMWF subset)
ax3 = axes[1, 0]
im3 = ax3.pcolormesh(
    ecmwf_subset['x'],
    ecmwf_subset['y'],
    ecmwf_subset['V10'],
    shading='auto',
    cmap='RdBu_r'
)
ax3.plot(lon_corners, lat_corners, 'k-', linewidth=2, label='SAR footprint')
ax3.set_xlabel('Longitude (°)')
ax3.set_ylabel('Latitude (°)')
ax3.set_title('ECMWF V10 (Before interpolation)')
ax3.legend()
plt.colorbar(im3, ax=ax3, label='V10 (m/s)')

# V10: After (interpolated on SAR grid)
ax4 = axes[1, 1]
im4 = ax4.scatter(
    sar_lon.flatten(),
    sar_lat.flatten(),
    c=result['V10'].values.flatten(),
    s=5,
    cmap='RdBu_r',
    vmin=im3.get_clim()[0],
    vmax=im3.get_clim()[1]
)
ax4.set_xlabel('Longitude (°)')
ax4.set_ylabel('Latitude (°)')
ax4.set_title('Interpolated V10 (After map_raster)')
plt.colorbar(im4, ax=ax4, label='V10 (m/s)')

plt.tight_layout()
plt.show()