# 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)

    [1]:

    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

    [2]:

    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]

    [3]:

    sar_dataset

    [3]:

    <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

    [4]:

    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]

    [5]:

    ecmwf_dataset

    [5]:

    <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

    [6]:

    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

    [7]:

    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)

    [8]:

    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)

    [9]:

    # 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)

    [10]:

    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

    [11]:

    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()

    [ ]: