Repository: SciTools/cartopy (python)

Issue:
NaturalEarthFeature is hard-coded to PlateCarree()
Currently the `Feature` class takes a `crs` argument at `__init__` time; `NaturalEarthFeature` hard-codes this as `PlateCarree()`.  This is generally fine, as it matches the Natural Earth data. It would be more flexible, though to pull this from the dataset if available, since shapefiles can optionally have a `.prj` file that contains the CRS in WKT format. We'd need to read that manually for pyshp, but Fiona makes this available as the `.crs` attribute on the base collection. I've confirmed that the Natural Earth data do include such `.prj` files.

This is relevant because downstream in MetPy we subclass `Feature` to provide US counties and states, (now) from the U.S. Census Bureau. These data are lat/lons defined on NAD83 (which uses GRS80) rather than WGS84. Not a huge difference, but always best to use what's in the data when available. The current API of `Feature` doesn't have a good way to set the correct feature at data load time.

Pull request description:
ENH: Use shapefile crs attribute when possible
When reading in the shapefile, use the .prj file information if it is present to set a crs attribute on the reader which can be used down the line by features to get the correct Projection.

Somewhat interesting that loading the wkt from the file produces a different CRS than the one going through Fiona's CRS first... I'm not sure if that is expected or not, so opening this PR as a draft for now.

closes ##906
closes #2289

Patch:
diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py
index 88c02856..d614afc4 100644
--- a/lib/cartopy/crs.py
+++ b/lib/cartopy/crs.py
@@ -696,6 +696,10 @@ class Projection(CRS, metaclass=ABCMeta):
             self.bounds = (x.min(), x.max(), y.min(), y.max())
             x0, x1, y0, y1 = self.bounds
             self.threshold = min(x1 - x0, y1 - y0) / 100.
+        elif self.is_geographic:
+            # If the projection is geographic without an area of use, assume
+            # the bounds are the full globe.
+            self.bounds = (-180, 180, -90, 90)
 
     @property
     def boundary(self):
diff --git a/lib/cartopy/feature/__init__.py b/lib/cartopy/feature/__init__.py
index 9804f5c7..f021b6b1 100644
--- a/lib/cartopy/feature/__init__.py
+++ b/lib/cartopy/feature/__init__.py
@@ -291,7 +291,10 @@ class NaturalEarthFeature(Feature):
             path = shapereader.natural_earth(resolution=self.scale,
                                              category=self.category,
                                              name=self.name)
-            geometries = tuple(shapereader.Reader(path).geometries())
+            reader = shapereader.Reader(path)
+            if reader.crs is not None:
+                self._crs = reader.crs
+            geometries = tuple(reader.geometries())
             _NATURAL_EARTH_GEOM_CACHE[key] = geometries
         else:
             geometries = _NATURAL_EARTH_GEOM_CACHE[key]
@@ -422,7 +425,10 @@ class GSHHSFeature(Feature):
                 # Load GSHHS geometries from appropriate shape file.
                 # TODO selective load based on bbox of each geom in file.
                 path = shapereader.gshhs(scale, level)
-                geoms = tuple(shapereader.Reader(path).geometries())
+                reader = shapereader.Reader(path)
+                if reader.crs is not None:
+                    self._crs = reader.crs
+                geoms = tuple(reader.geometries())
                 GSHHSFeature._geometries_cache[(scale, level)] = geoms
             for geom in geoms:
                 if extent is None or extent_geom.intersects(geom):
diff --git a/lib/cartopy/io/shapereader.py b/lib/cartopy/io/shapereader.py
index 704b2924..7312e3f5 100644
--- a/lib/cartopy/io/shapereader.py
+++ b/lib/cartopy/io/shapereader.py
@@ -32,10 +32,12 @@ import itertools
 from pathlib import Path
 from urllib.error import HTTPError
 
+from pyproj import CRS
 import shapefile
 import shapely.geometry as sgeom
 
 from cartopy import config
+import cartopy.crs as ccrs
 from cartopy.io import Downloader
 
 
@@ -136,6 +138,13 @@ class BasicReader:
     def __init__(self, filename, bbox=None, **kwargs):
         # Validate the filename/shapefile
         self._reader = reader = shapefile.Reader(filename, **kwargs)
+        # Try to set the CRS from the .prj file if it exists
+        self.crs = None
+        try:
+            with open(Path(filename).with_suffix('.prj'), 'rt') as fobj:
+                self.crs = ccrs.Projection(CRS.from_wkt(fobj.read()))
+        except FileNotFoundError:
+            pass
         self._bbox = bbox
         if reader.shp is None or reader.shx is None or reader.dbf is None:
             raise ValueError("Incomplete shapefile definition "
@@ -195,7 +204,13 @@ class FionaReader:
 
     def __init__(self, filename, bbox=None, **kwargs):
         self._data = []
-
+        # Try to set the CRS from the .prj file if it exists
+        self.crs = None
+        try:
+            with open(Path(filename).with_suffix('.prj'), 'rt') as fobj:
+                self.crs = ccrs.Projection(CRS.from_wkt(fobj.read()))
+        except FileNotFoundError:
+            pass
         with fiona.open(filename, **kwargs) as f:
             if bbox is not None:
                 assert len(bbox) == 4


Test patch:
diff --git a/lib/cartopy/tests/test_crs.py b/lib/cartopy/tests/test_crs.py
index bf433504..cffe2087 100644
--- a/lib/cartopy/tests/test_crs.py
+++ b/lib/cartopy/tests/test_crs.py
@@ -360,3 +360,25 @@ def test_transform_point_no_warning():
     with warnings.catch_warnings():
         warnings.simplefilter("error")
         p2.transform_point(1, 2, p)
+
+
+def test_geographic_bounds_no_area_of_use():
+    # Should default to global bounds if no area of use
+    wkt = ('GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",'
+           'SPHEROID["WGS_1984",6378137.0,298.257223563]],'
+           'PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]')
+    p = pyproj.CRS.from_wkt(wkt)
+    assert p.is_geographic
+    assert p.area_of_use is None
+    p_cartopy = ccrs.Projection(p)
+    assert p_cartopy.bounds == (-180, 180, -90, 90)
+
+
+def test_geographic_bounds_with_area_of_use():
+    # Should default to area of use bounds if available
+    p = pyproj.CRS.from_epsg(4267)
+    assert p.is_geographic
+    assert p.area_of_use is not None
+    p_cartopy = ccrs.Projection(p)
+    x0, x1, y0, y1 = p_cartopy.bounds
+    assert_array_equal((x1, y0, x0, y1), p.area_of_use.bounds)
diff --git a/lib/cartopy/tests/test_shapereader.py b/lib/cartopy/tests/test_shapereader.py
index 66c4e4d6..c13a3c75 100644
--- a/lib/cartopy/tests/test_shapereader.py
+++ b/lib/cartopy/tests/test_shapereader.py
@@ -62,6 +62,10 @@ class TestLakes:
         assert actual == expected
         assert lake_record.geometry == self.test_lake_geometry
 
+    def test_no_included_projection_file(self):
+        # No .prj file included with the lakes shapefile
+        assert self.reader.crs is None
+
     def test_bounds(self):
         if isinstance(self.reader, shp.BasicReader):
             # tests that a file which has a record with a bbox can
@@ -120,3 +124,15 @@ class TestRivers:
             if key in expected_attributes:
                 assert value == expected_attributes[key]
         assert river_record.geometry == self.test_river_geometry
+
+    def test_included_projection_file(self):
+        # This shapefile includes a .prj definition
+        wkt = ('GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",'
+            'ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],'
+            'ID["EPSG",6326]],PRIMEM["Greenwich",0,'
+            'ANGLEUNIT["Degree",0.0174532925199433]],CS[ellipsoidal,2],'
+            'AXIS["longitude",east,ORDER[1],'
+            'ANGLEUNIT["Degree",0.0174532925199433]],'
+            'AXIS["latitude",north,ORDER[2],'
+            'ANGLEUNIT["Degree",0.0174532925199433]]]')
+        assert self.reader.crs.to_wkt() == wkt