From Earth, a geostationary satellite looks like it is always in the same place, because it moves in the same direction and at the same rate the Earth spins. Image credit: NASA Solar System Exploration

Coordinate Types¶

Overview¶

Great circles use different types of coordinates when working with unit spheres and ellipsoids. This notebook will cover the different types of coordinates and how to convert between them.

Types of Coordinates
Convert Coordinates to All Coordinate Types
Plot Coordinates on a World Map

Prerequisites¶

Concepts	Importance	Notes
Numpy	Necessary	Used to work with large arrays
Pandas	Necessary	Used to read in and organize data (in particular dataframes)
Intro to Cartopy	Helpful	Will be used for adding maps to plotting
Matplotlib	Helpful	Will be used for plotting

Time to learn: 20 minutes

Imports¶

import numpy as np                                     # working with degrees and radians

import matplotlib.pyplot as plt                        # plotting a graph
from cartopy import crs as ccrs, feature as cfeature   # plotting a world map

Types of Coordinates¶

Geodesic Coordinates¶

Geodesic coordinates are latitude and longtiude and are measured from -90° South to 90° North and -180° East to 180° West measured from Greenwich.

Longitude lines are perpendicular to and latitude lines are parallel to the Equator from Wikipedia

Cartesian Coordinates¶

Cartesian coordinates describe points in space based on perpendicular axis lines that meet at a single point of origin, where any point’s position is described based on the distance to the origin along xyz axis.

A three dimensional Cartesian coordinate system, with origin O and axis lines X, Y and Z, oriented as shown by the arrows. The tick marks on the axes are one length unit apart. The black dot shows the point with coordinates x = 2, y = 3, and z = 4, or (2, 3, 4) from Wikipedia

Image Source: Three Dimensional Cartesian Coordinate System

Geodesic to Cartesian Coordinates

Assuming the Earth’s radius is 6378137 meters then:

x = radius * cos(latitude) * cos(longitude)

(1)

y = radius * cos(latitude) * sin(longitude)

(2)

z = radius * sin(latitude)

(3)

def cartesian_coordinates(latitude=None, longitude=None):
    earth_radius = 6378137  # meters
    latitude = np.deg2rad(latitude)
    longitude = np.deg2rad(longitude)
    cart_x = earth_radius * np.cos(latitude) * np.cos(longitude)
    cart_y = earth_radius * np.cos(latitude) * np.sin(longitude)
    cart_z = earth_radius * np.sin(latitude)
    return cart_x, cart_y, cart_z

Spherical Coordinates¶

Spherical coordinates describe points in space based on three values: radial distance (rho, r) along the radial line between point and the origin, polar angle (theta, θ) between the radial line and the polar axis, and azimuth angle (phi, φ) which is the angle of rotation of the radial line around the polar axis. With a fixed radius, the 3-point coordinates (r, θ, φ) provide a coordinate along a sphere.

Radial distance: distance from center to surface of sphere
Polar angle: angle between radial line and polar axis
Azimuth angle: angle around polar axis

Spherical Coordinate Description from Wikipedia

Image Source: Wikipedia - Spherical Coordinate System

Convert from cartesian (rectangular) coordinates spherical coordinates

ρ^2 = x^2 + y^2 + z^2

(4)

tan(θ) = \frac{y}{x}

(5)

φ = arccos(\frac{x}{\sqrt{x^2 + y^2 + z^2}})

(6)

Where, rho (ρ), theta (θ), phi (φ):

ρ = \sqrt{x^2 + y^2 + z^2}

(7)

θ = arctan(\frac{y}{x})

(8)

φ = arccos(\frac{x}{ρ})

(9)

def cartesian_to_spherical_coordinates(cart_x=None, cart_y=None, cart_z=None):
    rho = np.sqrt(cart_x**2 + cart_y**2 + cart_z**2)
    theta = np.arctan(cart_y/cart_x)
    phi = np.arccos(cart_z / rho)
    return rho, theta, phi

Polar Coordinates¶

Polar coordinates are a combination of latitude, longitude, and altitude from the center of the sphere (based on the radius).

Assuming the Earth’s radius is 6378137 meters then:

x = cos(latitude) * cos(longitude) * radius

(10)

y = cos(latitude) * sin(longitude) * radius

(11)

z = sin(latitude) * radius

(12)

def polar_coordinates(latitude=None, longitude=None):
    earth_radius = 6378137  # meters
    latitude = np.deg2rad(latitude)
    longitude = np.deg2rad(longitude)
    polar_x = np.cos(latitude) * np.sin(longitude) * earth_radius
    polar_y = np.cos(latitude) * np.cos(longitude) * earth_radius
    polar_z = np.sin(latitude) * earth_radius
    return polar_x, polar_y, polar_z

Convert City Coordinates to All Coordinate Types¶

Display Coordinates of Cities¶

First, we will read in the latitude and longitude coordinates from locations csv:

import pandas as pd

location_df = pd.read_csv("../location_coords.txt")
location_df = location_df.rename(columns=lambda x: x.strip()) # strip excess white space from column names and values
location_df

Add Columns for Additional Coordinate Types¶

location_df["cart_x"], location_df["cart_y"], location_df["cart_z"] = cartesian_coordinates(location_df["latitude"],
                                                                                            location_df["longitude"])
location_df["rho"], location_df["theta"], location_df["phi"] = cartesian_to_spherical_coordinates(location_df["cart_x"],
                                                                                                  location_df["cart_y"],
                                                                                                  location_df["cart_z"])
location_df["polar_x"], location_df["polar_y"], location_df["polar_z"] = polar_coordinates(location_df["latitude"],
                                                                                           location_df["longitude"])

location_df

# Save Output to a New Text File
location_df.to_csv("../location_full_coords.txt", index=False)

Plot Coordinates¶

World Map¶

Full world map from -180-180 and -90-90:

longitude east = 180

longitude west = -180

latitude north = 90

latitude south = -90

# Set up world map plot
fig = plt.subplots(figsize=(15, 10))
projection_map = ccrs.PlateCarree()
ax = plt.axes(projection=projection_map)
lon_west, lon_east, lat_south, lat_north = -180, 180, -90, 90
ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
ax.coastlines(color="black")
ax.add_feature(cfeature.BORDERS, edgecolor='grey')
ax.add_feature(cfeature.STATES, edgecolor="grey")
    
# Plot Latitude/Longitude Location
longitudes = location_df["longitude"] # longitude
latitudes = location_df["latitude"]   # latitude
plt.scatter(longitudes, latitudes, c="red")

plt.title("World Map with Locations")
plt.show()

/home/runner/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_1_states_provinces_lakes.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

---------------------------------------------------------------------------
error                                     Traceback (most recent call last)
File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/IPython/core/formatters.py:402, in BaseFormatter.__call__(self, obj)
    400     pass
    401 else:
--> 402     return printer(obj)
    403 # Finally look for special method names
    404 method = get_real_method(obj, self.print_method)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/IPython/core/pylabtools.py:170, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    167     from matplotlib.backend_bases import FigureCanvasBase
    168     FigureCanvasBase(fig)
--> 170 fig.canvas.print_figure(bytes_io, **kw)
    171 data = bytes_io.getvalue()
    172 if fmt == 'svg':

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/backend_bases.py:2155, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2152     # we do this instead of `self.figure.draw_without_rendering`
   2153     # so that we can inject the orientation
   2154     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2155         self.figure.draw(renderer)
   2156 if bbox_inches:
   2157     if bbox_inches == "tight":

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     92 @wraps(draw)
     93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94     result = draw(artist, renderer, *args, **kwargs)
     95     if renderer._rasterizing:
     96         renderer.stop_rasterizing()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
   3254             # ValueError can occur when resizing a window.
   3256     self.patch.draw(renderer)
-> 3257     mimage._draw_list_compositing_images(
   3258         renderer, self, artists, self.suppressComposite)
   3260     renderer.close_group('figure')
   3261 finally:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/mpl/geoaxes.py:524, in GeoAxes.draw(self, renderer, **kwargs)
    519         self.imshow(img, extent=extent, origin=origin,
    520                     transform=factory.crs, *factory_args[1:],
    521                     **factory_kwargs)
    522 self._done_img_factory = True
--> 524 return super().draw(renderer=renderer, **kwargs)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/axes/_base.py:3216, in _AxesBase.draw(self, renderer)
   3213 if artists_rasterized:
   3214     _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3216 mimage._draw_list_compositing_images(
   3217     renderer, self, artists, self.get_figure(root=True).suppressComposite)
   3219 renderer.close_group('axes')
   3220 self.stale = False

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    132 if not_composite or not has_images:
    133     for a in artists:
--> 134         a.draw(renderer)
    135 else:
    136     # Composite any adjacent images together
    137     image_group = []

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     68     if artist.get_agg_filter() is not None:
     69         renderer.start_filter()
---> 71     return draw(artist, renderer)
     72 finally:
     73     if artist.get_agg_filter() is not None:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/mpl/feature_artist.py:185, in FeatureArtist.draw(self, renderer)
    180     geoms = self._feature.geometries()
    181 else:
    182     # For efficiency on local maps with high resolution features (e.g
    183     # from Natural Earth), only create paths for geometries that are
    184     # in view.
--> 185     geoms = self._feature.intersecting_geometries(extent)
    187 stylised_paths = {}
    188 # Make an empty placeholder style dictionary for when styler is not
    189 # used.  Freeze it so that we can use it as a dict key.  We will need
    190 # to unfreeze all style dicts with dict(frozen) before passing to mpl.

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:309, in NaturalEarthFeature.intersecting_geometries(self, extent)
    302 """
    303 Returns an iterator of shapely geometries that intersect with
    304 the given extent.
    305 The extent is assumed to be in the CRS of the feature.
    306 If extent is None, the method returns all geometries for this dataset.
    307 """
    308 self.scaler.scale_from_extent(extent)
--> 309 return super().intersecting_geometries(extent)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:112, in Feature.intersecting_geometries(self, extent)
    109 if extent is not None and not np.isnan(extent[0]):
    110     extent_geom = sgeom.box(extent[0], extent[2],
    111                             extent[1], extent[3])
--> 112     return (geom for geom in self.geometries() if
    113             geom is not None and extent_geom.intersects(geom))
    114 else:
    115     return self.geometries()

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/feature/__init__.py:294, in NaturalEarthFeature.geometries(self)
    290 if key not in _NATURAL_EARTH_GEOM_CACHE:
    291     path = shapereader.natural_earth(resolution=self.scale,
    292                                      category=self.category,
    293                                      name=self.name)
--> 294     geometries = tuple(shapereader.Reader(path).geometries())
    295     _NATURAL_EARTH_GEOM_CACHE[key] = geometries
    296 else:

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/cartopy/io/shapereader.py:164, in BasicReader.geometries(self)
    152 def geometries(self):
    153     """
    154     Return an iterator of shapely geometries from the shapefile.
    155 
   (...)    162 
    163     """
--> 164     for shape in self._reader.iterShapes(bbox=self._bbox):
    165         # Skip the shape that can not be represented as geometry.
    166         if shape.shapeType != shapefile.NULL:
    167             yield sgeom.shape(shape)

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1483, in Reader.iterShapes(self, bbox)
   1479 if self.numShapes:
   1480     # Iterate exactly the number of shapes from shx header
   1481     for i in xrange(self.numShapes):
   1482         # MAYBE: check if more left of file or exit early? 
-> 1483         shape = self.__shape(oid=i, bbox=bbox)
   1484         if shape:
   1485             yield shape

File ~/micromamba/envs/cookbook-gc/lib/python3.13/site-packages/shapefile.py:1340, in Reader.__shape(self, oid, bbox)
   1338 # Read points - produces a list of [x,y] values
   1339 if nPoints:
-> 1340     flat = unpack("<%sd" % (2 * nPoints), f.read(16*nPoints))
   1341     record.points = list(izip(*(iter(flat),) * 2))
   1342 # Read z extremes and values

error: unpack requires a buffer of 336 bytes

<Figure size 1500x1000 with 2 Axes>

United States Map¶

Map of the United States roughly from -130 to -60 and 20 to 60:

longitude east = -60

longitude west = -130

latitude north = 60

latitude south = 20

# Set up United States map plot
fig = plt.subplots(figsize=(15, 10))
projection_map = ccrs.PlateCarree()
ax = plt.axes(projection=projection_map)
lon_west, lon_east, lat_south, lat_north = -130, -60, 20, 60
ax.set_extent([lon_west, lon_east, lat_south, lat_north], crs=projection_map)
ax.coastlines(color="black")
ax.add_feature(cfeature.BORDERS, edgecolor='grey')
ax.add_feature(cfeature.STATES, edgecolor="grey")
    
# Plot Latitude/Longitude Location
longitudes = location_df["longitude"] # longitude
latitudes = location_df["latitude"]   # latitude
plt.scatter(longitudes, latitudes, c="red")

plt.title("United States Map with Locations")
plt.show()

Summary¶

Coordinates on the Earth are measured in many different types of coordinate systems: Geodesic (latitude/longitude), cartesian, spherical, and polar. These coordinates will make future calculations simpler by converting a 2D coordinate like latitude/longitude into a 3D space that can be used for vector calculations.

In Python, coordinates can be mapped on to a world map via matplotlib and cartopy.

What’s next?¶

Great Circle arcs and paths

Resources and references¶

Foundations

Great Circle Terminology

Tutorials

Great Circle Arcs and Path