diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py
index 4e66588bb..4e975b4e9 100644
--- a/lib/cartopy/crs.py
+++ b/lib/cartopy/crs.py
@@ -644,27 +644,52 @@ class Projection(CRS, metaclass=ABCMeta):
     # Whether or not this projection can handle wrapped coordinates
     _wrappable = False
 
-    def __init__(self, *args, **kwargs):
+    def __init__(self, *args, bounds=None, transform_bounds=False, **kwargs):
+        """
+        Parameters
+        ----------
+        bounds: list or tuple, optional
+            Four element projection bounds (xmin, xmax, ymin, ymax).
+            If not provided defaults to "area_of_use" of the Coordinate
+            Reference System if it exists. This can be used to define an area
+            of use when one doesn't exist or to define a subset of the full
+            projection space.
+        transform_bounds: bool, optional
+            Transform bounds provided as lon(x)/lat(y) degrees to projection
+            space. This is a convenience when bounds are only known in degrees
+            for a non-geographic CRS. Defaults to False.
+
+        Additional arguments are passed directly to :class:`cartopy.crs.CRS`.
+
+        """
         super().__init__(*args, **kwargs)
-        self.bounds = None
-        if self.area_of_use:
+        if bounds is None and self.area_of_use is not None:
+            bounds = (
+                self.area_of_use.west,
+                self.area_of_use.east,
+                self.area_of_use.south,
+                self.area_of_use.north,
+            )
+            transform_bounds = True
+
+        if bounds is not None and transform_bounds:
             # Convert lat/lon bounds to projected bounds.
             # Geographic area of the entire dataset referenced to WGS 84
             # NB. We can't use a polygon transform at this stage because
             # that relies on the existence of the map boundary... the very
             # thing we're trying to work out! ;-)
-            x0 = self.area_of_use.west
-            x1 = self.area_of_use.east
-            y0 = self.area_of_use.south
-            y1 = self.area_of_use.north
+            x0, x1, y0, y1 = bounds
             lons = np.array([x0, x0, x1, x1])
             lats = np.array([y0, y1, y1, y0])
             points = self.transform_points(self.as_geodetic(), lons, lats)
             x = points[:, 0]
             y = points[:, 1]
-            self.bounds = (x.min(), x.max(), y.min(), y.max())
-            x0, x1, y0, y1 = self.bounds
-            self.threshold = min(x1 - x0, y1 - y0) / 100.
+            bounds = (x.min(), x.max(), y.min(), y.max())
+
+        self.bounds = None
+        if bounds is not None:
+            x0, x1, y0, y1 = self.bounds = tuple(bounds)
+            self.threshold = min(abs(x1 - x0), abs(y1 - y0)) / 100.
 
     @property
     def boundary(self):
diff --git a/lib/cartopy/tests/test_crs.py b/lib/cartopy/tests/test_crs.py
index 2d01015db..284d857db 100644
--- a/lib/cartopy/tests/test_crs.py
+++ b/lib/cartopy/tests/test_crs.py
@@ -332,3 +332,38 @@ def test_projection__from_string():
 
 def test_crs__from_pyproj_crs():
     assert ccrs.CRS(pyproj.CRS("EPSG:4326")) == "EPSG:4326"
+
+
+@pytest.mark.parametrize(
+    ("crs_input", "bounds", "transform_bounds", "exp_bounds", "exp_threshold"),
+    [
+        ("EPSG:4326", None, False,
+         (-180.0, 180.0, -90.0, 90.0), 1.8),
+        ("EPSG:4326", [-100.0, -90.0, 35.0, 45.0], False,
+         (-100.0, -90.0, 35.0, 45.0), 0.1),
+        ("EPSG:4326", [-100.0, -90.0, 35.0, 45.0], True,
+         (-100.0, -90.0, 35.0, 45.0), 0.1),
+        ("EPSG:6932", [-40000., 40000., -40000., 40000.], False,
+         (-40000., 40000., -40000., 40000.), 800.0),
+        ("EPSG:6932", [20.0, 25.0, -85.0, -78.0], True,
+         (190942.4609561831, 565329.6096711607,
+          505972.06807129766, 1257011.612161974), 3743.871487149776),
+        # Meteosat Second Generation (MSG) - SEVIRI - Flipped Geostationary
+        ({'proj': 'geos', 'lon_0': 0.0, 'a': 6378169.00, 'b': 6356583.80,
+          'h': 35785831.00, 'units': 'm'},
+         [5500000, -5500000, -5500000, 5500000], False,
+         (5500000, -5500000, -5500000, 5500000), 110000.0),
+        ({'proj': 'geos', 'lon_0': 0.0, 'a': 6378169.00, 'b': 6356583.80,
+          'h': 35785831.00, 'units': 'm'},
+         [-10.0, 10.0, -10.0, 10.0], True,
+         (-1084697.7494802547, 1084697.7494802547,
+          -1093480.6233566382, 1093480.6233566382), 21693.954989605096),
+    ]
+)
+def test_projection_with_bounds(crs_input, bounds, transform_bounds,
+                                exp_bounds, exp_threshold):
+    pcrs = pyproj.CRS.from_user_input(crs_input)
+    crs = ccrs.Projection(pcrs, bounds=bounds,
+                          transform_bounds=transform_bounds)
+    assert crs.bounds == exp_bounds
+    assert crs.threshold == exp_threshold