Possible vector rotation issue

[gauteh]

Code from Lenny

import numpy as np
import pyproj
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# -------------------------------------------------------
# Projection definitions
# -------------------------------------------------------

proj_from = pyproj.Proj("+proj=stere +lat_0=90 +lon_0=0 +lat_ts=70 +datum=WGS84")
proj_to = pyproj.Proj("+proj=latlong +datum=WGS84")

proj_from_crs = ccrs.Stereographic(
    central_latitude=90,
    central_longitude=0,
    true_scale_latitude=70
)

transformer = pyproj.Transformer.from_proj(proj_from, proj_to)
geod = pyproj.Geod(ellps="WGS84")

# -------------------------------------------------------
# Create a small 3x3 polar stereographic grid
# -------------------------------------------------------

x0 = 1000000 # 1000 km 
y0 = - 1000000 # 10000 km 
dx = 100000  # 100 km spacing

x = np.array([x0-dx, x0, x0+dx])
y = np.array([y0-dx, y0, y0+dx])

X, Y = np.meshgrid(x, y)

# -------------------------------------------------------
# Define grid unit vectors
# -------------------------------------------------------

# grid x-direction vector
u_x = np.ones_like(X)
v_x = np.zeros_like(X)

# grid y-direction vector
u_y = np.zeros_like(X)
v_y = np.ones_like(X)

# -------------------------------------------------------
# Rotation function (simplified version of repo code)
# -------------------------------------------------------

def rotate_vectors(reader_x, reader_y, u_component, v_component, proj_from, proj_to):

    delta_y = 10  # 10 meters northward in projection space

    transformer = pyproj.Transformer.from_proj(proj_from, proj_to)

    # transform base points
    lon, lat = transformer.transform(reader_x, reader_y)

    # transform northward shifted points
    lon2, lat2 = transformer.transform(reader_x, reader_y + delta_y)

    # compute azimuth of projected north
    azimuth = geod.inv(lon, lat, lon2, lat2)[0]

    rot_angle_vectors_rad = np.radians(azimuth)

    rot_angle_rad = -rot_angle_vectors_rad

    u_rot = (
        u_component * np.cos(rot_angle_rad)
        - v_component * np.sin(rot_angle_rad)
    )

    v_rot = (
        u_component * np.sin(rot_angle_rad)
        + v_component * np.cos(rot_angle_rad)
    )

    return u_rot, v_rot, lon, lat


# -------------------------------------------------------
# Rotate vectors
# -------------------------------------------------------

u_x_geo, v_x_geo, lon, lat = rotate_vectors(X, Y, u_x, v_x, proj_from, proj_to)
u_y_geo, v_y_geo, _, _ = rotate_vectors(X, Y, u_y, v_y, proj_from, proj_to)

# -------------------------------------------------------
# Plot grid vectors in stereographic coordinates
# -------------------------------------------------------

fig = plt.figure(figsize=(12,5))

pad = 50000 # 50 km
pad_deg = 2 # degrees

ax1 = fig.add_subplot(121,projection=proj_from_crs)
ax1.set_extent([x[0]-pad,x[-1]+pad,y[0]-pad ,y[-1]+pad],crs=proj_from_crs)
ax1.quiver(X, Y, u_x, v_x, color="blue", scale=10)
ax1.quiver(X, Y, u_y, v_y, color="red", scale=10)

ax1.set_title("Grid vectors in stereographic coordinates")
ax1.set_aspect("equal")

# -------------------------------------------------------
# Plot rotated vectors on lon/lat map
# -------------------------------------------------------

ax2 = fig.add_subplot(122, projection=ccrs.PlateCarree())

ax2.coastlines()
ax2.set_global()

ax2.set_extent([np.min(lon)-pad_deg,np.max(lon)+pad_deg,np.min(lat)-pad_deg ,np.max(lat)+pad_deg],crs=ccrs.PlateCarree())
# ax2.set_extent([x[0]-pad,x[-1]+pad,y[0]-pad ,y[-1]+pad],crs=proj_from_crs)

for seg in [(X[0,:],Y[0,:]), (X[-1,:],Y[-1,:]),
            (X[:,0],Y[:,0]), (X[:,-1],Y[:,-1])]:
    ax2.plot(seg[0], seg[1], color='blue', lw=1, transform=proj_from_crs)

ax2.quiver(
    lon, lat,
    u_x_geo, v_x_geo,
    color="blue",
    transform=ccrs.PlateCarree(),
    scale=10
)

ax2.quiver(
    lon, lat,
    u_y_geo, v_y_geo,
    color="red",
    transform=ccrs.PlateCarree(),
    scale=10
)
ax2.quiver(
    X, Y,
    u_x, v_x,
    color="black",
    transform=proj_from_crs,
    scale=10
)

ax2.quiver(
    X, Y,
    u_y, v_y,
    color="black",
    transform=proj_from_crs,
    scale=10
)

ax2.set_title("Vectors rotated into geographic coordinates")

plt.show()

[gauteh]

[lenny-frno]

`import numpy as np
import matplotlib.pyplot as plt
import pyproj
import cartopy.crs as ccrs
import cartopy.feature as cfeature

def run_projection_test(proj4, cartopy_crs, title):

# ── 1. SMALL 3×3 GRID ───────────────────────────────────────────────
x0 = 1_000_000
y0 = -1_000_000
dx = 100_000

x = np.array([x0-dx, x0, x0+dx])
y = np.array([y0-dx, y0, y0+dx])

X, Y = np.meshgrid(x, y)

dx = x[1] - x[0]
dy = y[1] - y[0]

# ── 2. PROJECTIONS ─────────────────────────────────────────────────
proj_crs = pyproj.CRS.from_proj4(proj4)
geo_crs  = pyproj.CRS.from_epsg(4326)

transformer = pyproj.Transformer.from_crs(proj_crs, geo_crs, always_xy=True)

lon_deg, lat_deg = transformer.transform(X, Y)

lon_rad = np.deg2rad(lon_deg)
lat_rad = np.deg2rad(lat_deg)

# ── 3. JACOBIAN ───────────────────────────────────────────────────
T = np.zeros((*X.shape, 2, 2))

T[...,0,0] = np.gradient(lon_rad, dx, axis=1)
T[...,0,1] = np.gradient(lon_rad, dy, axis=0)
T[...,1,0] = np.gradient(lat_rad, dx, axis=1)
T[...,1,1] = np.gradient(lat_rad, dy, axis=0)

# ── 4. METRIC MATRIX ─────────────────────────────────────────────
R = 6371008.8

M = np.zeros_like(T)
M[...,0,0] = R*np.cos(lat_rad)
M[...,1,1] = R

# ── 5. GRID BASIS VECTORS ───────────────────────────────────────
e_x = np.stack([np.ones_like(X), np.zeros_like(X)], axis=-1)
e_y = np.stack([np.zeros_like(X), np.ones_like(X)], axis=-1)

# ── 6. TRANSFORM VECTORS ───────────────────────────────────────
e_x_lonlat = np.rad2deg(np.einsum('...ij,...j->...i', T, e_x))
e_y_lonlat = np.rad2deg(np.einsum('...ij,...j->...i', T, e_y))

e_x_lonlat /= np.linalg.norm(e_x_lonlat, axis=-1, keepdims=True)
e_y_lonlat /= np.linalg.norm(e_y_lonlat, axis=-1, keepdims=True)

# ── 7. SANITY CHECK ─────────────────────────────────────────────
geo_phys_x = np.einsum('...ij,...j->...i', M, np.einsum('...ij,...j->...i', T, e_x))
geo_phys_y = np.einsum('...ij,...j->...i', M, np.einsum('...ij,...j->...i', T, e_y))

print(title)
print("‖M·J·e_x‖ ≈", np.linalg.norm(geo_phys_x,axis=-1))
print("‖M·J·e_y‖ ≈", np.linalg.norm(geo_phys_y,axis=-1))
print()

pc = ccrs.PlateCarree()

# ── 8. PLOT ─────────────────────────────────────────────────────
fig = plt.figure(figsize=(12,6))

pad = 100000

# Projection space
ax1 = fig.add_subplot(121, projection=cartopy_crs)

ax1.set_extent([x[0]-pad, x[-1]+pad, y[0]-pad, y[-1]+pad], crs=cartopy_crs)

ax1.add_feature(cfeature.LAND, facecolor='lightgray')
ax1.coastlines()
# grid outline
for seg in [(X[0,:],Y[0,:]), (X[-1,:],Y[-1,:]),
            (X[:,0],Y[:,0]), (X[:,-1],Y[:,-1])]:
    ax1.plot(seg[0], seg[1], transform=cartopy_crs, color="blue")
ax1.quiver(X, Y, e_x[...,0], e_x[...,1], transform=cartopy_crs, color="black", scale=10)
ax1.quiver(X, Y, e_y[...,0], e_y[...,1], transform=cartopy_crs, color="black", scale=10)

ax1.set_title(f"{title}\nProjection coordinates")

# Geographic space
ax2 = fig.add_subplot(122, projection=pc)

ax2.set_extent([x[0]-pad, x[-1]+pad, y[0]-pad, y[-1]+pad], crs=cartopy_crs)

ax2.add_feature(cfeature.LAND, facecolor='lightgray')
ax2.coastlines()
# grid outline
for seg in [(X[0,:],Y[0,:]), (X[-1,:],Y[-1,:]),
            (X[:,0],Y[:,0]), (X[:,-1],Y[:,-1])]:
    ax2.plot(seg[0], seg[1], transform=cartopy_crs, color="blue")
ax2.quiver(X, Y, e_x[...,0], e_x[...,1], transform=cartopy_crs, color="black", scale=10)
ax2.quiver(X, Y, e_y[...,0], e_y[...,1], transform=cartopy_crs, color="black", scale=10)

ax2.quiver(lon_deg, lat_deg, e_x_lonlat[...,0], e_x_lonlat[...,1],
           transform=pc, color="red", scale=10)

ax2.quiver(lon_deg, lat_deg, e_y_lonlat[...,0], e_y_lonlat[...,1],
           transform=pc, color="red", scale=10)

ax2.set_title("Vectors expressed in lon/lat")

ax2.gridlines(draw_labels=True)

plt.tight_layout()
plt.show()

─────────────────────────────────────────────────────────

CASE 1 : POLAR STEREOGRAPHIC

─────────────────────────────────────────────────────────

run_projection_test(
"+proj=stere +lat_0=90 +lon_0=-30 +lat_ts=90 +datum=WGS84",
ccrs.Stereographic(central_longitude=-30, central_latitude=90, true_scale_latitude=90),
"Polar stereographic"
)

─────────────────────────────────────────────────────────

CASE 2 : LAMBERT CONFORMAL

─────────────────────────────────────────────────────────

run_projection_test(
"+proj=lcc +lat_1=60 +lat_2=70 +lat_0=65 +lon_0=10 +datum=WGS84",
ccrs.LambertConformal(central_longitude=10, central_latitude=65,
standard_parallels=(60,70)),
"Lambert conformal"
)`

[lenny-frno]

[knutfrode]

Interesting, but what is the actual issue here?

…

On Fri, Mar 13, 2026, 17:33 lenny-frno ***@***.***> wrote: *lenny-frno* left a comment (OpenDrift/opendrift#1778) <#1778 (comment)> example.tranformation.vectors_LCC.png (view on web) <https://github.com/user-attachments/assets/b9889103-a377-49c2-8816-60d613eefc89> — Reply to this email directly, view it on GitHub <#1778 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABH25I6VTLB2VCXRYNPGQT34QQZ4BAVCNFSM6AAAAACWRGQM6GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHM2DANJWGQYTANRUGM> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[lenny-frno]

@knutfrode The code pasted by Gaute is a unit test of the transformation of the vectors (e_x,e_y) from the grid into (e_lon,e_lat) .
I just reproduced the logic implemented in opendrift where you compute locally the rotation to apply to each vector. The plot shows that in some cases , the simple rotation as implement in opendrift does not express correctly your vectors in terms of e_lon and e_lat. In the plot below Gaute's post the red and blue arrows should overlap the black arrows.
I suggested below, see my previous post, a more general solution based on the jacobian if the transformation xy --> lonlat. In the example of the application I'm taking the same projections (stereo and lambert) that were generating problems earlier. We see the overlap with some small error at the grid border because of the way the Jacobian is computed.

[knutfrode]

But is this not a plotting issue?
Unit vectors in the polarstereographic grid will point towards the next grid point when plotted in the native grid.
But when projected to lonlat, the same lines are no longer straight, and we do not anymore expect this alignment.
And also the new vector units are in degrees of longitude and latitude, which are not equal - more than a factor of 2 at these latitudes.

[lenny-frno]

@knutfrode. I don't think it is a plotting issue. Since after transformation of the native vector into the lon lat space we expect these transformed vectors (using a PlateCarree transform argument) to look the same as the native vectors (using a transform = projection) . That's what these lines are checking:

ax2.quiver(X, Y, e_x[...,0], e_x[...,1], transform=cartopy_crs, color="black", scale=10)
ax2.quiver(X, Y, e_y[...,0], e_y[...,1], transform=cartopy_crs, color="black", scale=10)

ax2.quiver(lon_deg, lat_deg, e_x_lonlat[...,0], e_x_lonlat[...,1],
           transform=pc, color="red", scale=10)

ax2.quiver(lon_deg, lat_deg, e_y_lonlat[...,0], e_y_lonlat[...,1],
           transform=pc, color="red", scale=10)

I agree that when projected to lonlat the lines between grid points are no longer straight especially if spaced by a lot but we can assume locally that we have a linear transformation. This depends on the grid spacing and gives a better result than just a rotation. By applying just a rotation you do not account for the potential stretching of the grid and you end up with a strong difference between the black arrows and the red/Blue arrows (see Gaute's comment on top).
For the unit point, This is what the metric matrix is for :

R = 6371008.8

M = np.zeros_like(T)
M[...,0,0] = R*np.cos(lat_rad)
M[...,1,1] = R

For the plot this is not used because I wanted to check the direction only. The conservation of norm is verified here :

# ── 7. SANITY CHECK ─────────────────────────────────────────────
geo_phys_x = np.einsum('...ij,...j->...i', M, np.einsum('...ij,...j->...i', T, e_x))
geo_phys_y = np.einsum('...ij,...j->...i', M, np.einsum('...ij,...j->...i', T, e_y))

print(title)
print("‖M·J·e_x‖ ≈", np.linalg.norm(geo_phys_x,axis=-1))
print("‖M·J·e_y‖ ≈", np.linalg.norm(geo_phys_y,axis=-1))
print()

Maybe some work needs to be done. Ideally the goal is to transform the native grid unit vectors to the lon lat vectors and keeping it unit free with a norm of one so that when we transform the winds or currents from native grid to lon lat space, without introducing any unit besides that the final wind / current vector should represent m/s.