How to use the cartopy.crs.Orthographic function in Cartopy

Ex7, Average lat=67.2, lon=-6.9
    """
    positions = np.array([(90, 0),
                          (60, 10),
                          (50, -20),
                          ])
    lats, lons = positions.T
    nvecs = lat_lon2n_E(rad(lats), rad(lons))

    # Find the horizontal mean position:
    n_EM_E = unit(np.sum(nvecs, axis=1).reshape((3, 1)))
    lat, lon = n_E2lat_lon(n_EM_E)
    lat, lon = deg(lat), deg(lon)
    print('Ex7, Average lat={0:2.1f}, lon={1:2.1f}'.format(lat[0], lon[0]))

    ax = plt.figure().gca(projection=ccrs.Orthographic(int(lon), int(lat)))
    ax.add_feature(cpf.OCEAN, zorder=0)
    ax.add_feature(cpf.LAND, zorder=0, edgecolor='black')
    ax.add_feature(cpf.COASTLINE)
    ax.add_feature(cpf.BORDERS, linestyle=':')
    ax.add_feature(cpf.LAKES, alpha=0.5)
    ax.add_feature(cpf.RIVERS)

    ax.set_global()
    ax.gridlines()

    # vector_crs = ccrs.Geodetic()
    vector_crs = ccrs.PlateCarree()

    ax.scatter(lon, lat, linewidth=5, marker='o', color='r',
               transform=vector_crs)
    ax.scatter(lons, lats, linewidth=5, marker='o', color='k',

def build_projections():  # {{{
    projections = {}
    projections['cyl'] = cartopy.crs.PlateCarree()
    projections['merc'] = cartopy.crs.Mercator()
    projections['mill'] = cartopy.crs.Miller()
    projections['mill2'] = cartopy.crs.Miller(central_longitude=180.)
    projections['moll'] = cartopy.crs.Mollweide()
    projections['moll2'] = cartopy.crs.Mollweide(central_longitude=180.)
    projections['robin'] = cartopy.crs.Robinson()
    projections['robin2'] = cartopy.crs.Robinson(central_longitude=180.)

    projections['ortho'] = cartopy.crs.Orthographic(central_longitude=-100.,
                                                    central_latitude=45.)
    projections['northpole'] = cartopy.crs.Orthographic(central_longitude=0.,
                                                        central_latitude=90.)
    projections['southpole'] = cartopy.crs.Orthographic(central_longitude=0.,
                                                        central_latitude=-90.)
    projections['atlantic'] = cartopy.crs.Orthographic(central_longitude=0.,
                                                       central_latitude=0.)
    projections['pacific'] = cartopy.crs.Orthographic(central_longitude=180.,
                                                      central_latitude=0.)
    projections['americas'] = cartopy.crs.Orthographic(central_longitude=-90.,
                                                       central_latitude=0.)
    projections['asia'] = cartopy.crs.Orthographic(central_longitude=90.,
                                                   central_latitude=0.)

    return projections  # }}}

z_djf_mean = cdutil.averager(z_djf, axis='t')
z_djf = z_djf - z_djf_mean
z_djf.id = 'z'

# Create an EOF solver to do the EOF analysis. Square-root of cosine of
# latitude weights are applied before the computation of EOFs.
solver = Eof(z_djf, weights='coslat')

# Retrieve the leading EOF, expressed as the covariance between the leading PC
# time series and the input SLP anomalies at each grid point.
eof1 = solver.eofsAsCovariance(neofs=1)

# Plot the leading EOF expressed as covariance in the European/Atlantic domain.
clevs = np.linspace(-75, 75, 11)
lons, lats = eof1.getLongitude()[:], eof1.getLatitude()[:]
proj = ccrs.Orthographic(central_longitude=-20, central_latitude=60)
ax = plt.axes(projection=proj)
ax.set_global()
ax.coastlines()
ax.contourf(lons, lats, eof1(squeeze=True).data, levels=clevs,
            cmap=plt.cm.RdBu_r, transform=ccrs.PlateCarree())
plt.title('EOF1 expressed as covariance', fontsize=16)

plt.show()

import harmonica as hm

# Load the global gravity grid
data = hm.datasets.fetch_gravity_earth()
print(data)

# Calculate normal gravity using the WGS84 ellipsoid
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell)
# The disturbance is the observed minus normal gravity (calculated at the
# observation point)
disturbance = data.gravity - gamma

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=160))
pc = disturbance.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="seismic"
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
ax.set_title("Gravity of disturbance of the Earth")
ax.coastlines()
plt.tight_layout()
plt.show()

Ratio between the central heights of the transformed image and original
        image.
    centrallonglat_xy : pandas.DataFrame
        xy position of the central longitude and latitude.
    """

    # If user doesn't provide Moon globe properties.
    if not iglobe:
        iglobe = ccrs.Globe(semimajor_axis=arad*1000.,
                            semiminor_axis=arad*1000., ellipse=None)

    # Set up Geodetic (long/lat), Plate Carree (usually long/lat, but not when
    # globe != WGS84) and Orthographic projections.
    geoproj = ccrs.Geodetic(globe=iglobe)
    iproj = ccrs.PlateCarree(globe=iglobe)
    oproj = ccrs.Orthographic(central_longitude=np.mean(llbd[:2]),
                              central_latitude=np.mean(llbd[2:]),
                              globe=iglobe)

    # Create and transform coordinates of image corners and edge midpoints.
    # Due to Plate Carree and Orthographic's symmetries, max/min x/y values of
    # these 9 points represent extrema of the transformed image.
    xll = np.array([llbd[0], np.mean(llbd[:2]), llbd[1]])
    yll = np.array([llbd[2], np.mean(llbd[2:]), llbd[3]])
    xll, yll = np.meshgrid(xll, yll)
    xll = xll.ravel()
    yll = yll.ravel()

    # [:,:2] because we don't need elevation data.
    res = iproj.transform_points(x=xll, y=yll, src_crs=geoproj)[:, :2]
    iextent = [min(res[:, 0]), max(res[:, 0]), min(res[:, 1]), max(res[:, 1])]

# Iris emits a warning due to the non-contiguous time dimension.
    warnings.simplefilter('ignore', UserWarning)
    z_djf_mean = z_djf.collapsed('time', iris.analysis.MEAN)
z_djf = z_djf - z_djf_mean

# Create an EOF solver to do the EOF analysis. Square-root of cosine of
# latitude weights are applied before the computation of EOFs.
solver = Eof(z_djf, weights='coslat')

# Retrieve the leading EOF, expressed as the covariance between the leading PC
# time series and the input SLP anomalies at each grid point.
eof1 = solver.eofsAsCovariance(neofs=1)

# Plot the leading EOF expressed as covariance in the European/Atlantic domain.
clevs = np.linspace(-75, 75, 11)
proj = ccrs.Orthographic(central_longitude=-20, central_latitude=60)
ax = plt.axes(projection=proj)
ax.coastlines()
ax.set_global()
iplt.contourf(eof1[0, 0], levels=clevs, cmap=plt.cm.RdBu_r)
ax.set_title('EOF1 expressed as covariance', fontsize=16)
plt.show()

This is the magnitude of the gravity vector of the Earth (gravitational
+ centrifugal) at 10 km height. The data is on a regular grid with 0.5 degree
spacing at 10km ellipsoidal height. It was generated from the spherical
harmonic model EIGEN-6C4 [Forste_etal2014]_.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import rockhound as rh

# Load the gravity grid
data = rh.fetch_gravity_earth()
print(data)

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=150))
pc = data.gravity.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Gravity of the Earth (EIGEN-6C4)")
ax.coastlines()
plt.tight_layout()
plt.show()

# Calculate normal gravity and the disturbance
ellipsoid = bl.WGS84
gamma = ellipsoid.normal_gravity(data.latitude, data.height_over_ell)
disturbance = data.gravity - gamma

# Reference the topography to the ellipsoid
topography_ell = data.topography + data.geoid

# Calculate the Bouguer planar correction and the topography-free disturbance.
# Use the default densities for the crust and ocean water.
bouguer = hm.bouguer_correction(topography_ell)
disturbance_topofree = disturbance - bouguer

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
pc = disturbance_topofree.plot.pcolormesh(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
)
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
ax.set_title("Topography-free (Bouguer) gravity of disturbance of the Earth")
ax.coastlines()
plt.tight_layout()
plt.show()

this case). Negative values indicate that the geoid is below the ellipsoid
surface and positive values that it is above. The data are on a regular grid
with 0.5 degree spacing and was generated from the spherical harmonic model
EIGEN-6C4 [Forste_etal2014]_.
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm

# Load the geoid grid
data = hm.datasets.fetch_geoid_earth()
print(data)

# Make a plot of data using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=100))
pc = data.geoid.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)
plt.colorbar(
    pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
)
ax.set_title("Geoid heights (EIGEN-6C4)")
ax.coastlines()
plt.tight_layout()
plt.show()

"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import harmonica as hm
import verde as vd

# Create a global computation grid with 1 degree spacing
region = [0, 360, -90, 90]
longitude, latitude = vd.grid_coordinates(region, spacing=1)

# Compute normal gravity for the WGS84 ellipsoid (the default) on its surface
gamma = hm.normal_gravity(latitude=latitude, height=0)

# Make a plot of the normal gravity using Cartopy
plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Orthographic())
ax.set_title("Normal gravity of the WGS84 ellipsoid")
ax.coastlines()
pc = ax.pcolormesh(longitude, latitude, gamma, transform=ccrs.PlateCarree())
plt.colorbar(
    pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.5
)
plt.tight_layout()
plt.show()