Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Short-cut when transforming PlateCarree -> PlateCarree in mpl. #260

Merged
merged 5 commits into from

2 participants

@pelson
Owner

Speeds up one PlateCarree contour from 15s to 8s.

A selection of performance test user-code will be available in the near future (not this PR).

@pelson pelson closed this
@pelson pelson reopened this
lib/cartopy/mpl/geoaxes.py
((5 lines not shown))
x_limits = projection.x_limits
- bypass = x.min() >= x_limits[0] and x.max() <= x_limits[1]
+ y_limits = projection.y_limits
+ bypass = (x.min() >= x_limits[0] and x.max() <= x_limits[1]
+ and y.min() > y_limits[0] and y.max() < y_limits[1])
@pelson Owner
pelson added a note

This change caused the image difference.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rhattersley
Owner

:+1: for this change in principle, but there are some detail issues to resolve.

lib/cartopy/mpl/geoaxes.py
((15 lines not shown))
+
+ src = self.source_projection
+ target = self.target_projection
+ mod = np.diff(target.x_limits)[0]
+
+ verts = src_path.vertices
+ x_lim = verts[:, 0].min(), verts[:, 0].max()
+ y_lim = verts[:, 1].min(), verts[:, 1].max()
+
+ potential = (src.y_limits[0] <= y_lim[0] and
+ src.y_limits[1] >= y_lim[1])
+
+ try:
+ bboxes, proj_offset = target._bbox_and_offset(src)
+ except (ValueError, TypeError):
+ potential = False
@rhattersley Owner

The main issue I have is with the overall structure of this section. The use of the try..except, the private method, and the projection-specific logic (some of which you inherited) all suggest the code is in the wrong place.

I think it would be better if the InterProjectionTransform class knew nothing of _CylindricalProjection or PlateCarree projections. Instead, it should enquire of its projection "Do I really need to transform these vertices? If not, should I use some simple alternative coordinate values instead?".

For example:

bypass_vertices = self.source_projection.can_bypass(self.target_projection, src_path.vertices)
if bypass_vertices is not None:
    if bypass_vertices is not src_path.vertices:
        src_path = mpath.Path(bypass_vertices, src_path.codes)
    return src_path
@pelson Owner
pelson added a note

I'd be happy enough with such an interface, and agree with your evaluation that the code is in the wrong place. Do you want me to do that in this PR?

@rhattersley Owner

Yes please. :smile:

@rhattersley Owner

... even if it's just for the PlateCarree subset.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pelson
Owner

@rhattersley - I've now updated this PR.

@rhattersley
Owner

I've now updated this PR.

That's much nicer - thank you :smile:

I'll have a proper look through now...

lib/cartopy/crs.py
((8 lines not shown))
+ to be transformed to wrap appropriately.
+
+ >>> import cartopy.crs as ccrs
+ >>> src = ccrs.PlateCarree(central_longitude=10)
+ >>> bboxes, offset = ccrs.PlateCarree()._bbox_and_offset(src)
+ >>> print bboxes
+ [[-180, -170.0], [-170.0, 180]]
+ >>> print offset
+ 10.0
+
+ The returned values are longitudes in ``other_plate_carree``'s
+ coordinate system.
+
+ """
+ if not isinstance(other_plate_carree, PlateCarree):
+ raise TypeError('pc_proj2 must be a PlateCaree instance.')
@rhattersley Owner

Since this is now a proper private method it'd be simpler to drop this check now. (But if it must stay, pc_proj2 looks to be a legacy term and should be other_plate_carree.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
lib/cartopy/crs.py
((19 lines not shown))
+ coordinate system.
+
+ """
+ if not isinstance(other_plate_carree, PlateCarree):
+ raise TypeError('pc_proj2 must be a PlateCaree instance.')
+
+ self_params = self.proj4_params.copy()
+ other_params = other_plate_carree.proj4_params.copy()
+ self_params.pop('lon_0'), other_params.pop('lon_0')
+ if self_params != other_params:
+ raise ValueError('All proj4 params (other than lon_0) of '
+ '"other_plate_carree" must be equal to '
+ 'those of self.')
+
+ central_lon_0 = self.proj4_params['lon_0']
+ central_lon_1 = other_plate_carree.proj4_params['lon_0']
@rhattersley Owner

There's a confusing double meaning for _0 here, especially when combined as lon_0. I suggest you stick with the self vs other convention you've used earlier. For example:

  • central_lon_0 -> self_lon_0
  • central_lon_0_offset -> lon_0_offset
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
lib/cartopy/crs.py
((58 lines not shown))
+
+ # Optimise the PlateCarree -> PlateCarree case where no
+ # wrapping or interpolation needs to take place.
+ if return_value is None and isinstance(src_crs, PlateCarree):
+ mod = np.diff(src_crs.x_limits)[0]
+
+ x_lim = vertices[:, 0].min(), vertices[:, 0].max()
+ y_lim = vertices[:, 1].min(), vertices[:, 1].max()
+
+ potential = (self.y_limits[0] <= y_lim[0] and
+ self.y_limits[1] >= y_lim[1])
+
+ try:
+ bboxes, proj_offset = self._bbox_and_offset(src_crs)
+ except (ValueError, TypeError):
+ potential = False
@rhattersley Owner

Using a try..except block here still seems odd when it's really just normal logic flow. Perhaps the structure in this function should be:

if PlateCarree -> PlateCarree and all params matching except lon_0:
    Compute bboxes and proj_offset (now with no failure modes)
    Check if vertices sit inside bboxes (as before)
@pelson Owner
pelson added a note

I'm not keen on this. If anybody else wants to make use of _bbox_and_offset they will also have to do this checking - also, without running the same tests in _bbox_and_offset there is no guarantee that its result is sensible. It is for that reason that I'm comfortable with a try...except block inline here.
Obviously, if you insist I will go ahead and implement the same checks in both methods to remove the try...except - but I'm not keen for _bbox_and_offset to avoid raising Exceptions.

@rhattersley Owner

The bit I really care about is you shouldn't catch the TypeError. You've already checked the type so you don't expect a TypeError to be raised, so if one does crop up it's a genuine error and shouldn't be silenced.

After that, I still think the use of the ValueError suggests a more elegant refactor is possible ... but I'm not going to get hung up on that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
lib/cartopy/tests/mpl/test_axes.py
((7 lines not shown))
+
+ # of the 3 paths, 2 of them cannot be short-cutted.
+ pth1 = mpath.Path([[0.5, 0], [10, 10]])
+ pth2 = mpath.Path([[0.5, 91], [10, 10]])
+ pth3 = mpath.Path([[-0.5, 0], [10, 10]])
+
+ trans = InterProjectionTransform(src, target)
+
+ counter = CallCounter(target, 'project_geometry')
+
+ with counter:
+ trans.transform_path(pth1)
+ trans.transform_path(pth2)
+ trans.transform_path(pth3)
+
+ assert_equal(counter.count, 2)
@rhattersley Owner

By grouping them we don't know which one used the short-cut.

@pelson Owner
pelson added a note

No, but it doesn't really matter. I don't even test that the result is correct. It is purely an "is this optimised" test. There are other tests which will pick up badly projected polygons.

@pelson Owner
pelson added a note

FWIW I've updated this test to be more specific.

@rhattersley Owner

Thanks - that clarifies which paths should be optimised, which might be handy if/when we ever come revisit this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
lib/cartopy/tests/test_crs.py
((10 lines not shown))
+ ([[-180, -170], [-170, 180]], 10),
+ ([[-180, 170], [170, 180]], -10),
+ ([[-180, 180], [180, 180]], 360),
+ ([[-180, -180], [-180, 180]], -360),
+ ]
+
+ assert len(target) == len(central_lons)
+
+ for expected, (s_lon0, t_lon0) in zip(target, central_lons):
+ expected_bboxes, expected_offset = expected
+
+ src = ccrs.PlateCarree(central_longitude=s_lon0)
+ target = ccrs.PlateCarree(central_longitude=t_lon0)
+
+ bbox, offset = src._bbox_and_offset(target)
+ print s_lon0, t_lon0, bbox, offset
@rhattersley Owner

cough

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pelson
Owner

I've done all actions bar the one I've commented on. Lets thrash this out and merge in the next 2 hours :wink:.

Thanks for the review @rhattersley.

@pelson
Owner

Just added a commit which addresses the try except.

@rhattersley
Owner

Nice work @pelson :smile:

It's faster and it's cleaned out some hacky code from the InterProjectionTransform. :+1:

@rhattersley rhattersley merged commit 767b85d into SciTools:master
@pelson
Owner

Thanks @rhattersley .

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
10 docs/source/whats_new.rst
@@ -1,3 +1,13 @@
+What's new in Cartopy 0.8
+*************************
+
+:Release: 0.8.0
+:Date: unreleased
+
+* Phil Elson added a major performance improvement when plotting data from PlateCarree onto a
+ PlateCarree map. (:pull:`260`)
+
+
What's new in Cartopy 0.7
*************************
View
106 lib/cartopy/crs.py
@@ -466,6 +466,31 @@ def _rings_to_multi_polygon(self, rings, is_ccw):
multi_poly = sgeom.MultiPolygon()
return multi_poly
+ def quick_vertices_transform(self, vertices, src_crs):
+ """
+ Where possible, return a vertices array transformed to this CRS from
+ the given vertices array of shape ``(n, 2)`` and the source CRS.
+
+ .. important::
+
+ This method may return None to indicate that the vertices cannot
+ be transformed quickly, and a more complex geometry transformation
+ is required (see :meth:`cartopy.crs.Projection.project_geometry`).
+
+ """
+ return_value = None
+
+ if self == src_crs:
+ x = vertices[:, 0]
+ y = vertices[:, 1]
+ x_limits = self.x_limits
+ y_limits = self.y_limits
+ if (x.min() >= x_limits[0] and x.max() <= x_limits[1]
+ and y.min() >= y_limits[0] and y.max() <= y_limits[1]):
+ return_value = vertices
+
+ return return_value
+
class _RectangularProjection(Projection):
"""
@@ -511,6 +536,87 @@ def __init__(self, central_longitude=0.0):
def threshold(self):
return 0.5
+ def _bbox_and_offset(self, other_plate_carree):
+ """
+ Returns a pair of (xmin, xmax) pairs and an offset which can be used
+ for identification of whether data in ``other_plate_carree`` needs
+ to be transformed to wrap appropriately.
+
+ >>> import cartopy.crs as ccrs
+ >>> src = ccrs.PlateCarree(central_longitude=10)
+ >>> bboxes, offset = ccrs.PlateCarree()._bbox_and_offset(src)
+ >>> print bboxes
+ [[-180, -170.0], [-170.0, 180]]
+ >>> print offset
+ 10.0
+
+ The returned values are longitudes in ``other_plate_carree``'s
+ coordinate system.
+
+ .. important::
+
+ The two CRSs must be identical in every way, other than their
+ central longitudes. No checking of this is done.
+
+ """
+ self_lon_0 = self.proj4_params['lon_0']
+ other_lon_0 = other_plate_carree.proj4_params['lon_0']
+
+ lon_0_offset = other_lon_0 - self_lon_0
+
+ lon_lower_bound_0 = self.x_limits[0]
+ lon_lower_bound_1 = (other_plate_carree.x_limits[0] + lon_0_offset)
+
+ if lon_lower_bound_1 < self.x_limits[0]:
+ lon_lower_bound_1 += np.diff(self.x_limits)[0]
+
+ lon_lower_bound_0, lon_lower_bound_1 = sorted(
+ [lon_lower_bound_0, lon_lower_bound_1])
+
+ bbox = [[lon_lower_bound_0, lon_lower_bound_1],
+ [lon_lower_bound_1, lon_lower_bound_0]]
+
+ bbox[1][1] += np.diff(self.x_limits)[0]
+
+ return bbox, lon_0_offset
+
+ def quick_vertices_transform(self, vertices, src_crs):
+ return_value = super(PlateCarree,
+ self).quick_vertices_transform(vertices, src_crs)
+
+ # Optimise the PlateCarree -> PlateCarree case where no
+ # wrapping or interpolation needs to take place.
+ if return_value is None and isinstance(src_crs, PlateCarree):
+ self_params = self.proj4_params.copy()
+ src_params = src_crs.proj4_params.copy()
+ self_params.pop('lon_0'), src_params.pop('lon_0')
+
+ xs, ys = vertices[:, 0], vertices[:, 1]
+
+ potential = (self_params == src_params and
+ self.y_limits[0] <= ys.min() and
+ self.y_limits[1] >= ys.max())
+
+ if potential:
+ mod = np.diff(src_crs.x_limits)[0]
+ bboxes, proj_offset = self._bbox_and_offset(src_crs)
+ x_lim = xs.min(), xs.max()
+ y_lim = ys.min(), ys.max()
+ for poly in bboxes:
+ # Arbitrarily choose the number of moduli to look
+ # above and below the -180->180 range. If data is beyond
+ # this range, we're not going to transform it quickly.
+ for i in [-1, 0, 1, 2]:
+ offset = mod * i - proj_offset
+ if ((poly[0] + offset) <= x_lim[0]
+ and (poly[1] + offset) >= x_lim[1]):
+ return_value = vertices + [[-offset, 0]]
+ break
+ if return_value is not None:
+ break
+
+ return return_value
+
class TransverseMercator(_RectangularProjection):
def __init__(self, central_longitude=0.0):
View
85 lib/cartopy/mpl/geoaxes.py
@@ -26,13 +26,12 @@
import matplotlib
import matplotlib.axes
-import matplotlib.collections as mcollections
from matplotlib.image import imread
import matplotlib.transforms as mtransforms
import matplotlib.patches as mpatches
import matplotlib.path as mpath
import numpy as np
-import shapely.geometry
+import shapely.geometry as sgeom
from cartopy import config
import cartopy.crs as ccrs
@@ -139,15 +138,15 @@ def transform_path_non_affine(self, src_path):
if result is not None:
return result
- bypass = self.source_projection == self.target_projection
- if bypass:
- projection = self.source_projection
- if isinstance(projection, ccrs._CylindricalProjection):
- x = src_path.vertices[:, 0]
- x_limits = projection.x_limits
- bypass = x.min() >= x_limits[0] and x.max() <= x_limits[1]
- if bypass:
- return src_path
+ # Allow the vertices to be quickly transformed, if
+ # quick_vertices_transform allows it.
+ new_vertices = self.target_projection.quick_vertices_transform(
+ src_path.vertices, self.source_projection)
+ if new_vertices is not None:
+ if new_vertices is src_path.vertices:
+ return src_path
+ else:
+ return mpath.Path(new_vertices, src_path.codes)
if src_path.vertices.shape == (1, 2):
return mpath.Path(self.transform(src_path.vertices))
@@ -440,9 +439,9 @@ def _get_extent_geom(self, crs=None):
x1, x2 = self.projection.x_limits
y1, y2 = self.projection.y_limits
- domain_in_src_proj = shapely.geometry.Polygon([[x1, y1], [x2, y1],
- [x2, y2], [x1, y2],
- [x1, y1]])
+ domain_in_src_proj = sgeom.Polygon([[x1, y1], [x2, y1],
+ [x2, y2], [x1, y2],
+ [x1, y1]])
# Determine target projection based on requested CRS.
if crs is None:
@@ -466,7 +465,7 @@ def _get_extent_geom(self, crs=None):
' coordinate system {!r}'.format(crs))
# Calculate intersection with boundary and project if necesary.
- boundary_poly = shapely.geometry.Polygon(self.projection.boundary)
+ boundary_poly = sgeom.Polygon(self.projection.boundary)
if proj != self.projection:
# Erode boundary by threshold to avoid transform issues.
# This is a workaround for numerical issues at the boundary.
@@ -493,9 +492,9 @@ def set_extent(self, extents, crs=None):
# plt.ylim - allowing users to set None for a minimum and/or
# maximum value
x1, x2, y1, y2 = extents
- domain_in_crs = shapely.geometry.LineString([[x1, y1], [x2, y1],
- [x2, y2], [x1, y2],
- [x1, y1]])
+ domain_in_crs = sgeom.LineString([[x1, y1], [x2, y1],
+ [x2, y2], [x1, y2],
+ [x1, y1]])
r = self.projection.project_geometry(domain_in_crs, crs)
x1, y1, x2, y2 = r.bounds
@@ -616,56 +615,6 @@ def set_yticks(self, ticks, minor=False, crs=None):
return super(GeoAxes, self).set_yticks(yticks, minor)
-# def geod_circle_meters(self, lon_0, lat_0, radius, npts=80, **kwargs):
-# # radius is in meters
-# geod = self.projection.as_geodetic()
-#
-# az = np.linspace(0, 360, npts)
-# lats = np.zeros(npts) + lat_0
-# lons = np.zeros(npts) + lon_0
-# distances = np.zeros(npts) + radius
-#
-# lons, lats, _reverse_az = geod.fwd(lons, lats, az, distances,
-# radians=False)
-# ll = np.concatenate([lons[:, None], lats[:, None]], 1)
-# from matplotlib.patches import Polygon
-# poly = Polygon(ll, transform=cartopy.prj.PlateCarree(), **kwargs)
-# self.add_patch(poly)
-# return poly
-#
-# def gshhs_line(self, outline_color='k', domain=None,
-# resolution='low', **kwargs):
-# # domain is a shapely geometry (Polygon or MultiPolygon)
-# import cartopy.gshhs as gshhs
-## import cartopy.spherical as spherical
-# from matplotlib.collections import PatchCollection, LineCollection
-#
-# paths = []
-#
-# projection = self.projection
-#
-# if domain is None:
-# domain = self.map_domain(ccrs.PlateCarree())
-#
-# for points in gshhs.read_gshhc(gshhs.fnames[resolution],
-# poly=False, domain=domain):
-# paths.extend(patch.geos_to_path(
-# shapely.geometry.LineString(points))
-# )
-#
-## slinestring = shapely.geometry.LineString(points)
-## projected = projection.project_geometry(slinestring)
-## paths.extend(patch.geos_to_path(projected))
-#
-# collection = PatchCollection([mpatches.PathPatch(pth)
-# for pth in paths],
-# edgecolor=outline_color, facecolor='none',
-# transform=ccrs.PlateCarree(),
-# **kwargs
-# )
-#
-# self.add_collection(collection, autolim=False)
-
def stock_img(self, name='ne_shaded'):
"""
Add a standard image to the map.
View
36 lib/cartopy/tests/mpl/test_axes.py
@@ -17,10 +17,16 @@
import unittest
+
+import matplotlib.path as mpath
import matplotlib.pyplot as plt
+from nose.tools import assert_equal
import numpy as np
+
import cartopy.crs as ccrs
+from cartopy.mpl.geoaxes import InterProjectionTransform
+from .test_caching import CallCounter
class TestNoSpherical(unittest.TestCase):
@@ -49,5 +55,33 @@ def test_pcolormesh(self):
self.ax.pcolormesh(self.data, transform=ccrs.Geodetic())
+def test_transform_PlateCarree_shortcut():
+ src = ccrs.PlateCarree(central_longitude=0)
+ target = ccrs.PlateCarree(central_longitude=180)
+
+ # of the 3 paths, 2 of them cannot be short-cutted.
+ pth1 = mpath.Path([[0.5, 0], [10, 10]])
+ pth2 = mpath.Path([[0.5, 91], [10, 10]])
+ pth3 = mpath.Path([[-0.5, 0], [10, 10]])
+
+ trans = InterProjectionTransform(src, target)
+
+ counter = CallCounter(target, 'project_geometry')
+
+ with counter:
+ trans.transform_path(pth1)
+ # pth1 should allow a short-cut.
+ assert_equal(counter.count, 0)
+
+ with counter:
+ trans.transform_path(pth2)
+ assert_equal(counter.count, 1)
+
+ with counter:
+ trans.transform_path(pth3)
+ assert_equal(counter.count, 2)
+
+
if __name__ == '__main__':
- unittest.main()
+ import nose
+ nose.runmodule(argv=['-s', '--with-doctest'], exit=False)
View
27 lib/cartopy/tests/test_crs.py
@@ -22,6 +22,7 @@
import numpy as np
from numpy.testing import assert_array_almost_equal as assert_arr_almost_eq
+from nose.tools import assert_equal
import cartopy.crs as ccrs
@@ -130,6 +131,32 @@ def test_pickle():
assert pc == ccrs.PlateCarree()
+def test_PlateCarree_shortcut():
+ central_lons = [[0, 0], [0, 180], [0, 10], [10, 0], [-180, 180], [
+ 180, -180]]
+
+ target = [([[-180, -180], [-180, 180]], 0),
+ ([[-180, 0], [0, 180]], 180),
+ ([[-180, -170], [-170, 180]], 10),
+ ([[-180, 170], [170, 180]], -10),
+ ([[-180, 180], [180, 180]], 360),
+ ([[-180, -180], [-180, 180]], -360),
+ ]
+
+ assert len(target) == len(central_lons)
+
+ for expected, (s_lon0, t_lon0) in zip(target, central_lons):
+ expected_bboxes, expected_offset = expected
+
+ src = ccrs.PlateCarree(central_longitude=s_lon0)
+ target = ccrs.PlateCarree(central_longitude=t_lon0)
+
+ bbox, offset = src._bbox_and_offset(target)
+
+ assert_equal(offset, expected_offset)
+ assert_equal(bbox, expected_bboxes)
+
+
if __name__ == '__main__':
import nose
nose.runmodule(argv=['-s', '--with-doctest'], exit=False)
Something went wrong with that request. Please try again.