reconstruct pole prototype

[larsbuntemeyer]

prototype code to reconstruct rotated pole if missing:

def reconstruct_pole(rlon, rlat, lon, lat):

    def sph2cart(lon, lat):
        lon = np.deg2rad(lon)
        lat = np.deg2rad(lat)

        x = np.cos(lat) * np.cos(lon)
        y = np.cos(lat) * np.sin(lon)
        z = np.sin(lat)

        return np.stack([x, y, z], axis=-1)
    
    RLON, RLAT = np.meshgrid(rlon, rlat)

    v_rot = sph2cart(RLON, RLAT).reshape(-1,3)
    v_geo = sph2cart(lon, lat).reshape(-1,3)

    H = v_rot.T @ v_geo
    U, S, Vt = np.linalg.svd(H)
    R = Vt.T @ U.T
    north_pole_vec = R @ np.array([0,0,1])

    lon_pole = np.degrees(np.arctan2(north_pole_vec[1], north_pole_vec[0]))
    lat_pole = np.degrees(np.arcsin(north_pole_vec[2]))

    return lon_pole, lat_pole

ds = cx.domain("EUR-12")
reconstruct_pole(ds.rlon, ds.rlat, ds.lon, ds.lat)