ENH: Use shapefile crs attribute when possible

lib/cartopy/feature/__init__.py

@@ -291,7 +291,10 @@ def geometries(self):
             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 @@ def intersecting_geometries(self, extent):
                 # 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):

lib/cartopy/io/shapereader.py

@@ -32,10 +32,12 @@
 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.CRS(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.CRS(CRS.from_wkt(fobj.read()))
+        except FileNotFoundError:
+            pass
         with fiona.open(filename, **kwargs) as f:
             if bbox is not None:
                 assert len(bbox) == 4

lib/cartopy/tests/test_shapereader.py

@@ -62,6 +62,10 @@ def test_record(self):
         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 @@ def test_record(self):
             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