Merge pull request #98 from djhoese/bugfix-longlat-crs

Remove special handling of geographic (longlat) CRSes

pycoast/cw_base.py

@@ -29,9 +29,9 @@
 from typing import Generator
 
 import numpy as np
-import pyproj
 import shapefile
 from PIL import Image
+from pyproj import CRS, Proj
 
 try:
     from pyresample import AreaDefinition
@@ -45,16 +45,11 @@ def get_resolution_from_area(area_def):
     """Get the best resolution for an area definition."""
     x_size = area_def.width
     y_size = area_def.height
-    prj = Proj(area_def.crs if hasattr(area_def, "crs") else area_def.proj_str)
-    if prj.is_latlong():
-        x_ll, y_ll = prj(area_def.area_extent[0], area_def.area_extent[1])
-        x_ur, y_ur = prj(area_def.area_extent[2], area_def.area_extent[3])
-        x_resolution = (x_ur - x_ll) / x_size
-        y_resolution = (y_ur - y_ll) / y_size
-    else:
-        x_resolution = (area_def.area_extent[2] - area_def.area_extent[0]) / x_size
-        y_resolution = (area_def.area_extent[3] - area_def.area_extent[1]) / y_size
+    x_resolution = abs(area_def.area_extent[2] - area_def.area_extent[0]) / x_size
+    y_resolution = abs(area_def.area_extent[3] - area_def.area_extent[1]) / y_size
     res = min(x_resolution, y_resolution)
+    if "degree" in area_def.crs.axis_info[0].unit_name:
+        res = _estimate_metered_resolution_from_degrees(area_def.crs, res)
 
     if res > 25000:
         return "c"
@@ -68,6 +63,12 @@ def get_resolution_from_area(area_def):
         return "f"
 
 
+def _estimate_metered_resolution_from_degrees(crs: CRS, resolution_degrees: float) -> float:
+    major_radius = crs.datum.ellipsoid.semi_major_metre
+    # estimate by taking the arc length using the radius
+    return major_radius * math.radians(resolution_degrees)
+
+
 class _CoordConverter:
     """Convert coordinates from one space to in-bound image pixel column and row.
 
@@ -91,7 +92,8 @@ def __init__(self, coord_ref: str, area_def: AreaDefinition):
             raise ValueError(f"'coord_ref' must be one of {pretty_coord_refs}.")
         self._convert_method = convert_methods[coord_ref]
 
-    def _check_area_def(self, area_def):
+    @staticmethod
+    def _check_area_def(area_def):
         if AreaDefinition is None:
             raise ImportError(
                 "Missing required 'pyresample' module, please "
@@ -130,16 +132,6 @@ def hash_dict(dict_to_hash: dict) -> str:
     return dhash.hexdigest()
 
 
-class Proj(pyproj.Proj):
-    """Wrapper around pyproj to add in 'is_latlong'."""
-
-    def is_latlong(self):
-        if hasattr(self, "crs"):
-            return self.crs.is_geographic
-        # pyproj<2.0
-        return super(Proj, self).is_latlong()
-
-
 class ContourWriterBase(object):
     """Base class for contourwriters. Do not instantiate.
 
@@ -1058,16 +1050,10 @@ def _get_bounding_box_lonlat_sides(area_extent, x_size, y_size, prj):
     x_range = np.linspace(x_ll, x_ur, num=x_size)
     y_range = np.linspace(y_ll, y_ur, num=y_size)
 
-    if prj.is_latlong():
-        lons_s1, lats_s1 = x_ll * np.ones(y_range.size), y_range
-        lons_s2, lats_s2 = x_range, y_ur * np.ones(x_range.size)
-        lons_s3, lats_s3 = x_ur * np.ones(y_range.size), y_range
-        lons_s4, lats_s4 = x_range, y_ll * np.ones(x_range.size)
-    else:
-        lons_s1, lats_s1 = prj(np.ones(y_range.size) * x_ll, y_range, inverse=True)
-        lons_s2, lats_s2 = prj(x_range, np.ones(x_range.size) * y_ur, inverse=True)
-        lons_s3, lats_s3 = prj(np.ones(y_range.size) * x_ur, y_range, inverse=True)
-        lons_s4, lats_s4 = prj(x_range, np.ones(x_range.size) * y_ll, inverse=True)
+    lons_s1, lats_s1 = prj(np.ones(y_range.size) * x_ll, y_range, inverse=True)
+    lons_s2, lats_s2 = prj(x_range, np.ones(x_range.size) * y_ur, inverse=True)
+    lons_s3, lats_s3 = prj(np.ones(y_range.size) * x_ur, y_range, inverse=True)
+    lons_s4, lats_s4 = prj(x_range, np.ones(x_range.size) * y_ll, inverse=True)
     return (lons_s1, lons_s2, lons_s3, lons_s4), (lats_s1, lats_s2, lats_s3, lats_s4)
 
 
@@ -1092,12 +1078,7 @@ def _get_pixel_index(shape, area_extent, x_size, y_size, prj, x_offset=0, y_offs
     shape_data = np.array(shape.points if hasattr(shape, "points") else shape)
     lons = shape_data[:, 0]
     lats = shape_data[:, 1]
-
-    if prj.is_latlong():
-        x_ll, y_ll = prj(area_extent[0], area_extent[1])
-        x_ur, y_ur = prj(area_extent[2], area_extent[3])
-    else:
-        x_ll, y_ll, x_ur, y_ur = area_extent
+    x_ll, y_ll, x_ur, y_ur = area_extent
 
     x, y = prj(lons, lats)
 

pycoast/tests/test_pycoast.py

@@ -28,6 +28,7 @@
 import pytest
 import shapefile
 from PIL import Image, ImageFont
+from pyproj import CRS
 from pyresample.geometry import AreaDefinition
 from pytest_lazyfixture import lazy_fixture
 
@@ -153,6 +154,20 @@ def dateline_2():
 DATELINE_2 = dateline_2()
 
 
+def lonlat_0(pm=0):
+    """Create longlat projection over Cuba."""
+    proj4_string = "+proj=longlat +lon_0=0.0 +ellps=WGS84"
+    if pm:
+        proj4_string += f" +pm={pm}"
+    area_extent = (pm + -90.0, 15.0, pm + -70.0, 25.0)
+    area_def = (proj4_string, area_extent)
+    return area_def
+
+
+LONLAT_0 = lonlat_0()
+LONLAT_0_180 = lonlat_0(pm=180)
+
+
 def nh():
     """Create the nh area."""
     proj4_string = "+proj=laea +lat_0=90 +lon_0=0 +a=6371228.0 +units=m"
@@ -444,18 +459,45 @@ def test_geos(self, cw_pil, new_test_image):
                     lat_placement="lr",
                 ),
             ),
+            (
+                "lonlat_boundary_cross.png",
+                (640, 480),
+                LONLAT_0,
+                4,
+                dict(
+                    fill="blue",
+                    write_text=False,
+                    outline="blue",
+                    minor_outline="blue",
+                    lon_placement="b",
+                    lat_placement="lr",
+                ),
+            ),
+            (
+                "lonlat_boundary_cross.png",
+                (640, 480),
+                LONLAT_0_180,
+                4,
+                dict(
+                    fill="blue",
+                    write_text=False,
+                    outline="blue",
+                    minor_outline="blue",
+                    lon_placement="b",
+                    lat_placement="lr",
+                ),
+            ),
         ],
     )
     @pytest.mark.usefixtures("cd_test_dir")
     def test_grid(self, cw_pil, new_test_image, filename, shape, area_def, level, grid_kwargs):
-        grid_img = Image.open(os.path.join(LOCAL_DIR, filename))
-
         img = new_test_image("RGB", shape, filename)
 
         cw_pil.add_coastlines(img, area_def, resolution="l", level=level)
         font = pil_font_16
         cw_pil.add_grid(img, area_def, (10.0, 10.0), (2.0, 2.0), font=font, **grid_kwargs)
 
+        grid_img = Image.open(os.path.join(LOCAL_DIR, filename))
         assert images_match(grid_img, img), "Writing of grid failed"
 
     def test_grid_nh(self, cw_pil, new_test_image):
@@ -1770,12 +1812,43 @@ def test_add_grid_from_dict_agg(self):
         res = np.array(img)
         assert fft_metric(grid_data, res), "Writing grid from dict agg failed"
 
+    def test_lonlat_pm_change(self):
+        """Test that a longlat projection with a non-0 prime meridian is handled correctly."""
+        from pyresample.geometry import AreaDefinition
+
+        from pycoast import ContourWriterAGG
+
+        area_def1 = AreaDefinition("", "", "", "+proj=longlat", 640, 480, (-55.0, -35.0, -5.0, 3.0))
+        area_def2 = AreaDefinition("", "", "", "+proj=longlat +pm=180", 640, 480, (-55.0, -35.0, -5.0, 3.0))
+        area_def3 = AreaDefinition("", "", "", "+proj=longlat +pm=180", 640, 480, (125.0, -35.0, 175.0, 3.0))
+        img1 = Image.new("RGBA", (640, 480))
+        img2 = Image.new("RGBA", (640, 480))
+        img3 = Image.new("RGBA", (640, 480))
+
+        cw = ContourWriterAGG(gshhs_root_dir)
+        cw.add_coastlines(img1, area_def1, resolution="l", level=4)
+        cw.add_coastlines(img2, area_def2, resolution="l", level=4)
+        cw.add_coastlines(img3, area_def3, resolution="l", level=4)
+
+        # with a prime meridian shift and extents shift, the images should be the same
+        np.testing.assert_allclose(np.array(img1), np.array(img3))
+        # with only a prime meridian shift and same extents, the images should be completely different
+        _assert_all_notclose(np.array(img1), np.array(img2))
+        # with only an extents change, the images should be completely different
+        _assert_all_notclose(np.array(img2), np.array(img3))
+
+
+def _assert_all_notclose(*args, **kwargs):
+    with pytest.raises(AssertionError):
+        np.testing.assert_allclose(*args, **kwargs)
+
 
 class FakeAreaDef:
     """A fake area definition object."""
 
     def __init__(self, proj4_string, area_extent, x_size, y_size):
-        self.proj_str = self.proj_dict = self.crs = proj4_string
+        self.proj_str = self.proj_dict = proj4_string
+        self.crs = CRS.from_user_input(proj4_string)
         self.area_extent = area_extent
         self.width = x_size
         self.height = y_size
@@ -2003,16 +2076,32 @@ def test_cache_nocache_consistency(
         # no cache version and manual version should be the same
         np.testing.assert_allclose(np.array(background_img3), np.array(background_img4), atol=0)
 
-    def test_get_resolution(self):
+    @pytest.mark.parametrize(
+        ("crs_input", "extents", "size", "exp_resolution"),
+        [
+            (
+                "+proj=stere +lon_0=8.00 +lat_0=50.00 +lat_ts=50.00 +ellps=WGS84",
+                (-3363403.31, -2291879.85, 2630596.69, 2203620.1),
+                (640, 480),
+                "l",
+            ),
+            (
+                "+proj=stere +lon_0=8.00 +lat_0=50.00 +lat_ts=50.00 +ellps=WGS84",
+                (-3363403.31, -2291879.85, 2630596.69, 2203620.1),
+                (6400, 4800),
+                "h",
+            ),
+            ("+proj=longlat +ellps=WGS84", (-45.0, -35.0, 5.0, 3.0), (640, 480), "l"),
+            ("+proj=longlat +pm=180 +ellps=WGS84", (-45.0, -35.0, 5.0, 3.0), (640, 480), "l"),
+        ],
+    )
+    def test_get_resolution(self, crs_input, extents, size, exp_resolution):
         """Get the automagical resolution computation."""
         from pycoast import get_resolution_from_area
 
-        proj4_string = "+proj=stere +lon_0=8.00 +lat_0=50.00 +lat_ts=50.00 +ellps=WGS84"
-        area_extent = (-3363403.31, -2291879.85, 2630596.69, 2203620.1)
-        area_def = FakeAreaDef(proj4_string, area_extent, 640, 480)
-        assert get_resolution_from_area(area_def) == "l"
-        area_def = FakeAreaDef(proj4_string, area_extent, 6400, 4800)
-        assert get_resolution_from_area(area_def) == "h"
+        crs = CRS.from_user_input(crs_input)
+        area_def = FakeAreaDef(crs, extents, *size)
+        assert get_resolution_from_area(area_def) == exp_resolution
 
 
 def _create_polygon_shapefile(fn: pathlib.Path, polygon_coords: list) -> None: