From 8784761df5827e768c57d4eca6b367fbad2e843e Mon Sep 17 00:00:00 2001 From: Elliott Sales de Andrade Date: Sat, 13 Oct 2018 04:07:07 -0400 Subject: [PATCH] Fix projected line strings getting truncated early. The limit to prevent an infinite loop attempting to split line segments "small enough" was applied to the entire line string, which meant the entire path could get truncated. Instead, make the limit apply to each segment individually. Since each segment is somewhat smaller, drop the limit a bit to compensate for having to go through more segments. This keeps things running in a reasonable time. --- lib/cartopy/_trace.cpp | 5 ++-- lib/cartopy/tests/test_polygon.py | 45 +++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+), 2 deletions(-) diff --git a/lib/cartopy/_trace.cpp b/lib/cartopy/_trace.cpp index c8d2f4c95..e13de3947 100644 --- a/lib/cartopy/_trace.cpp +++ b/lib/cartopy/_trace.cpp @@ -1,5 +1,5 @@ /* -# (C) British Crown Copyright 2010 - 2016, Met Office +# (C) British Crown Copyright 2010 - 2018, Met Office # # This file is part of cartopy. # @@ -563,7 +563,8 @@ void _project_segment(GEOSContextHandle_t handle, t_current = 0.0; state = get_state(p_current, gp_domain, handle); - while(t_current < 1.0 && lines.size() < 500) + size_t old_lines_size = lines.size(); + while(t_current < 1.0 && (lines.size() - old_lines_size) < 100) { //std::cerr << "Bisecting" << std::endl; #ifdef DEBUG diff --git a/lib/cartopy/tests/test_polygon.py b/lib/cartopy/tests/test_polygon.py index 7d4086d02..8d7900024 100644 --- a/lib/cartopy/tests/test_polygon.py +++ b/lib/cartopy/tests/test_polygon.py @@ -18,6 +18,7 @@ from __future__ import (absolute_import, division, print_function) import numpy as np +import pytest import shapely.geometry as sgeom import shapely.wkt @@ -135,6 +136,50 @@ def test_project_previous_infinite_loop(self): src = ccrs.PlateCarree() src._attach_lines_to_boundary(multi_line_strings, True) + @pytest.mark.parametrize('proj', + [ccrs.InterruptedGoodeHomolosine, ccrs.Mollweide]) + def test_infinite_loop_bounds(self, proj): + # test a polygon which used to get stuck in an infinite loop but is now + # erroneously clipped. + # see https://github.com/SciTools/cartopy/issues/1131 + + # These names are for IGH; effectively the same for Mollweide. + bottom = [0., 70.] + right = [0., 90.] + top = [-180., 90.] + left = [-180., 70.] + verts = np.array([ + bottom, + right, + top, + left, + bottom, + ]) + bad_path = sgeom.Polygon(verts) + + target = proj() + source = ccrs.PlateCarree() + + projected = target.project_geometry(bad_path, source) + + # When transforming segments was broken, the resulting path did not + # close, and either filled most of the domain, or a smaller portion + # than it should. Check that the bounds match the individual points at + # the expected edges. + projected_left = target.transform_point(left[0], left[1], source) + assert projected.bounds[0] == pytest.approx(projected_left[0], + rel=target.threshold) + projected_bottom = target.transform_point(bottom[0], bottom[1], source) + assert projected.bounds[1] == pytest.approx(projected_bottom[1], + rel=target.threshold) + projected_right = target.transform_point(right[0], right[1], source) + assert projected.bounds[2] == pytest.approx(projected_right[0], + rel=target.threshold, + abs=1e-8) + projected_top = target.transform_point(top[0], top[1], source) + assert projected.bounds[3] == pytest.approx(projected_top[1], + rel=target.threshold) + def test_3pt_poly(self): projection = ccrs.OSGB() polygon = sgeom.Polygon([(-1000, -1000),