2

I'm currently trying to produce maps of a shoreline, which runs roughly from south-east - north-west in an UTM projection to visualize shoreline change. By default Cartopy maps are in north up. Is there any way to rotate the map in a way that the coastline is parallel to the x-axis? This would save me a lot of space, since I have to create lots of subplots.

Some code to play with:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt

subplot_extent = [557565, 562329, 6070049, 6073791]

request = cimgt.OSM()
projection = ccrs.UTM(32)

fig = plt.figure(figsize=(9,7))
ax = fig.add_subplot(projection=projection)
ax.set_extent(subplot_extent, projection)
ax.add_image(request, 13)

This produces this map with north up: Map with North up

Expected output: Map with coast parallel to x axis

1
  • This is a bit old, but I really want to know this also!
    – Ajean
    Oct 30, 2023 at 20:37

1 Answer 1

2

I can't say how to do it with a UTM projection, but I was able to do something similar using a custom-modified Oblique Mercator projection where I've enabled setting the "+gamma" parameter to proj4 (normally if not specified it's set to match alpha so that north is always up).

The class looks almost exactly like the ObliqueMercator class (as of Cartopy 0.22), with that extra gamma specification. It appears to work nicely, at least in the conditions I've tested so far. Intuitively setting the right azimuth and gamma turns out to be a bit tricky for a particular desired rotation, but in this case just setting gamma to 0 works.

Doing exactly what you did above but using the custom projection with a central lat/lon and rotation of 35 degrees:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt

# Either put the definition of the custom class here or import it from another file

request = cimgt.OSM()
projection = RotatedObliqueMercator(central_latitude=54.790877, central_longitude=9.932319, azimuth=35, gamma=0)

fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
ax.set_extent((-3000, 3000, -1200, 1200), crs=projection)
ax.add_image(request, 13)
gdl = ax.gridlines(draw_labels=True)

Rotated Output I used bounds of (-3000, 3000, -1200, 1200) in the custom projection coordinates, plus added a gridliner just to see if it worked.  

Here's the code for the class I used:

from cartopy.crs import Projection, Mercator
import shapely.geometry as sgeom

class RotatedObliqueMercator(Projection):
    """An Oblique Mercator projection that maintains rotation."""

    def __init__(self, central_longitude=0.0, central_latitude=0.0,
                 false_easting=0.0, false_northing=0.0,
                     scale_factor=1.0, azimuth=0.0, globe=None, gamma=0.0):
        """
        Parameters
        ----------
        central_longitude: optional
            The true longitude of the central meridian in degrees.
            Defaults to 0.
        central_latitude: optional
            The true latitude of the planar origin in degrees. Defaults to 0.
        false_easting: optional
            X offset from the planar origin in metres. Defaults to 0.
        false_northing: optional
            Y offset from the planar origin in metres. Defaults to 0.
        scale_factor: optional
            Scale factor at the central meridian. Defaults to 1.
        azimuth: optional
            Azimuth of centerline clockwise from north at the center point of
            the centre line. Defaults to 0.
        globe: optional
            An instance of :class:`cartopy.crs.Globe`. If omitted, a default
            globe is created.

        Notes
        -----
        The 'Rotated Mercator' projection can be achieved using Oblique
        Mercator with `azimuth` ``=90``.
        """
        if np.isclose(azimuth, 90):
            # Exactly 90 causes coastline 'folding'.
            azimuth -= 1e-3

        # The only difference between this and the built-in ObliqueMercator is that we have
        # hard-coded gamma to be 0 so as to maintain the alpha rotation
        proj4_params = [('proj', 'omerc'), ('lonc', central_longitude),
                        ('lat_0', central_latitude), ('k', scale_factor),
                        ('x_0', false_easting), ('y_0', false_northing),
                        ('alpha', azimuth), ('units', 'm'), ('gamma', gamma)]

        super().__init__(proj4_params, globe=globe)

        # Couple limits to those of Mercator - delivers acceptable plots, and
        #  Mercator has been through much more scrutiny.
        mercator = Mercator(
            central_longitude=central_longitude,
            globe=globe,
            false_easting=false_easting,
            false_northing=false_northing,
            scale_factor=scale_factor,
        )
        self._x_limits = mercator.x_limits
        self._y_limits = mercator.y_limits
        self.threshold = mercator.threshold

    @property
    def boundary(self):
        x0, x1 = self.x_limits
        y0, y1 = self.y_limits
        return sgeom.LinearRing([(x0, y0), (x0, y1),
                                 (x1, y1), (x1, y0),
                                 (x0, y0)])

    @property
    def x_limits(self):
        return self._x_limits

    @property
    def y_limits(self):
        return self._y_limits

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.