Merge pull request #80 from CardiacModelling/sea-fixes-and-checks

Sea fixes and checks

README.md

@@ -138,6 +138,13 @@ A a great map without BNG coordinates can be found at https://opentopomap.org.
 
 The sea is a bit messy in these files, as the values depend on mean sea level in each 10x10 km^2 area (OS Tile) relative to OS datum (0m) level [which is mean sea level in Newlyn, Cornwall](https://en.wikipedia.org/wiki/Ordnance_datum).
 
+### The "sea"
+
+During preprocessing, the `nevis` module makes an attempt to classify points as "sea" or "not sea".
+This is done by assessing which below sea-level points can be accessed from the edge of the map, without ever going above sea-level.
+As a result, the mouths of small rivers, and large sections of wide rivers, are classified as "sea".
+All "sea" points are assigned height data to create a gentle slope pointing back towards the land.
+
 ### Hill tops
 
 Names of hill and mountain tops are taken from [The Database of British and Irish Hills v17.2](http://www.hills-database.co.uk), which is made available under a CC-BY license.

interpolation.py

@@ -20,16 +20,12 @@
 r = nevis.spacing()
 r2 = r // 2             # Height is taken at center of grid square
 
-# Get interpolants
-name = 'interpolation'
-f = nevis.linear_interpolant()
-g = nevis.spline(verbose=True)
-
 #
 # Squares and objects to draw on the maps
 #
-# Squares as (square, lines), where each square is (xlo, xhi, ylo, yhi) and
-# each line is (x0, x1, y0, y1)
+# Squares as (square, lines, sea), where each square is (xlo, xhi, ylo, yhi),
+# where each line is (x0, x1, y0, y1), and where sea is a bool indicating
+# whether the sea mask should be shown
 squares = []
 # Labels as name : Coords
 labels = {}
@@ -56,7 +52,7 @@
 lines.append([214675, 216675, 771275, 771275])
 # Diagonal, from valley up to peak
 lines.append([216675 - 1400, 216675, 771275 - 1400, 771275])
-squares.append((square, lines))
+squares.append((square, lines, False))
 
 #
 # The (firth of) clyde near Dumbarton Rock NS 4001 7391
@@ -69,7 +65,7 @@
 lines.append([239925, 239925, 674475, 672525])
 # Near Rock, horizontal
 lines.append([239925, 239225, 674475, 674475])
-squares.append((square, lines))
+squares.append((square, lines, True))
 
 #
 # Margate, near Dreamland and its Scenic Railway coaster 635107, 170534
@@ -81,11 +77,7 @@
 lines.append([635125, 635125, 170225, 171825])
 # Through a river and into the sea
 lines.append([633125, 636125, 161825, 161825])
-squares.append((square, lines))
-
-#
-# Teignmouth beach SX 9441 7301
-#
+squares.append((square, lines, True))
 
 #
 # Sea issue NT68: Coast near North Berwick 359353, 682985
@@ -96,10 +88,10 @@
 lines.append([360025, 369975, 684025, 684025])
 # Through square, vertically
 lines.append([365025, 365025, 680025, 689975])
-squares.append((square, lines))
+squares.append((square, lines, True))
 # Through below-sea level "land" in square below
 lines.append([362825, 362825, 677925, 680525])
-squares.append((square, lines))
+squares.append((square, lines, True))
 
 #
 # Sea issues NR35,24,34,44,33: Islay, south-east coast
@@ -114,7 +106,7 @@
 lines.append([129025, 129025, 649975, 647025])
 lines.append([120025, 149975, 645025, 645025])
 lines.append([130025, 139975, 635025, 635025])
-squares.append((square, lines))
+squares.append((square, lines, True))
 
 #
 # Sea issues: NR56,57: Jura south-east coast
@@ -126,7 +118,7 @@
 lines.append([156025, 159975, 672025, 672025])
 lines.append([153025, 159975, 666025, 666025])
 lines.append([150025, 159975, 661525, 661525])
-squares.append((square, lines))
+squares.append((square, lines, True))
 
 #
 # Sea issues: NR65,64,74,76,75: North-west Kintyre and south-west Knapdale
@@ -139,7 +131,7 @@
 lines.append([160025, 175025, 655025, 655025])
 lines.append([160025, 169975, 645025, 645025])
 lines.append([165025, 165025, 659975, 640025])
-squares.append((square, lines))
+squares.append((square, lines, True))
 
 #
 # End of the world
@@ -149,24 +141,51 @@
 lines.append([-1025, 1025, 5525, 5525])
 lines.append([5525, 5525, -1025, 1025])
 lines.append([-1025, 1025, -1025, 1025])
-squares.append((square, lines))
+squares.append((square, lines, False))
 
 w, h = nx * r, ny * r
 square = [w - 8000, w + 8000, -8000, 8000]
 lines = []
 lines.append([w - 1025, w + 1025, 5525, 5525])
 lines.append([w - 5525, w - 5525, -1025, 1025])
 lines.append([w - 1025, w + 1025, 1025, -1025])
-squares.append((square, lines))
+squares.append((square, lines, False))
 
 w, h = nx * r, ny * r
 square = [w - 8000, w + 8000, h - 8000, h + 8000]
 lines = []
 lines.append([w - 1025, w + 1025, h - 5525, h - 5525])
 lines.append([w - 5525, w - 5525, h - 1025, h + 1025])
 lines.append([w - 1025, w + 1025, h - 1025, h + 1025])
-squares.append((square, lines))
+squares.append((square, lines, False))
 
+#
+# Teignmouth
+#
+square = [290000, 300000, 65000, 80000]
+lines = []
+# Horizontally through squares
+lines.append([293025, 295025, 72525, 72525])
+squares.append((square, lines, True))
+
+#
+# Holme fen 520483, 289083
+#
+square = [510000, 570000, 280000, 320000]
+lines = []
+# Horizontally through squares
+lines.append([520025, 531025, 289075, 289075])
+lines.append([558825, 560825, 309825, 309825])
+lines.append([559825, 559825, 308825, 310825])
+squares.append((square, lines, True))
+
+#
+# Get interpolants
+#
+name = 'interpolation'
+f = nevis.linear_interpolant()
+g = nevis.spline(verbose=True)
+h = lambda x, y: 1 if nevis.is_sea(nevis.Coords(x, y)) else 0
 
 #
 # Ensure results directory exists
@@ -183,7 +202,7 @@
     zoom=0.03,
     headless=True)
 for ii, sq in enumerate(squares):
-    square, line = sq
+    square, line, show_sea_mask = sq
     x0, x1, y0, y1 = square
     x0, y0 = ff(x0, y0)
     x1, y1 = ff(x1, y1)
@@ -200,7 +219,7 @@
 # Figure 2: Zoomed maps
 #
 for ii, sq in enumerate(squares):
-    square, lines = sq
+    square, lines, show_sea_mask = sq
     fig, ax, data, ff = nevis.plot(
         boundaries=square,
         labels=labels,
@@ -226,11 +245,13 @@
     fig.savefig(path)
 
     # Figure 3: Line from known top to you
+    fs = (f, g, h) if show_sea_mask else (f, g)
+
     for jj, line in enumerate(lines):
         x0, x1, y0, y1 = line
         p0, p1 = nevis.Coords(x0, y0), nevis.Coords(x1, y1)
         fig, ax, q0, q1 = nevis.plot_line(
-            (f, g), p0, p1, figsize=(14, 10), headless=True)
+            fs, p0, p1, figsize=(14, 10), evaluations=1000, headless=True)
         ax.minorticks_on()
         ax.grid(which='major')
         ax.grid(which='minor', color='#eeeeee')

nevis/_bng.py

@@ -489,11 +489,12 @@ def pub(name=None):
 Coords.ben = Coords(gridx=216666, gridy=771288)
 Coords.fen = Coords(gridx=520483, gridy=289083)
 Coords.pub = {
-    #'Ye olde trip to jerusalem': Coords(gridx=457034, gridy=339443},
     'Bear': Coords(gridx=451473, gridy=206135),
     'Canal house': Coords(gridx=457307, gridy=339326),
     'MacSorleys': Coords(gridx=258809, gridy=665079),
     'Sheffield tap': Coords(gridx=435847, gridy=387030),
+    'Ship Inn': Coords(gridx=293945, gridy=72712),
+    'Laurel Inn': Coords(gridx=495227, gridy=505037),
 }
 
 

nevis/_nevis_version.py

@@ -1,5 +1,5 @@
 #
 # Version info
 #
-__version_tuple__ = (0, 0, 8)
+__version_tuple__ = (0, 1, 0)
 __version__ = '.'.join(str(x) for x in __version_tuple__)

nevis/_os_terrain_50.py

@@ -136,8 +136,12 @@ def is_sea(coords):
     if _heights is None:
         gb()
     x, y = coords.grid // nevis.spacing()
-    return _heights[y, x] < 0 and (
-        (y * _heights.shape[1] + x) not in _not_sea)
+    try:
+        if _heights[y, x] >= 0:
+            return False
+    except IndexError:  # Anything off-map counts as sea
+        return True
+    return (y * _heights.shape[1] + x) not in _not_sea
 
 
 def download_os_terrain_50(force=False):

nevis/_plot.py

@@ -289,10 +289,6 @@ def plot_line(f, point_1, point_2, label_1='Point 1', label_2='Point 2',
     """
     Draws a line between two points and evaluates a function along it.
 
-    Note that this method assumes you will want to write the figure to disk
-    with :meth:`fig.savefig()`. If you want to display it using ``pyplot``,
-    set ``headless=False``.
-
     Arguments:
 
     ``f``
@@ -312,7 +308,8 @@ def plot_line(f, point_1, point_2, label_1='Point 1', label_2='Point 2',
     ``figsize``
         The default figure size
     ``headless``
-        Set to ``True`` to create the figure without using pyplot.
+        Set to ``True`` to create the figure without using pyplot, i.e. to save
+        with ``figsave`` instead of calling ``plt.show()``.
     ``verbose``
         Set to ``True`` to print updates to ``stdout`` while plotting.
 
@@ -339,10 +336,11 @@ def plot_line(f, point_1, point_2, label_1='Point 1', label_2='Point 2',
         fs = [f]
     else:
         fs = f
-        for f in fs:
+        for i, f in enumerate(fs):
             if not callable(f):
                 raise ValueError(
-                    'f must be a callable or a sequence of callables.')
+                    'f must be a callable or a sequence of callables: found'
+                    f' non-callable at index {i}.')
 
     # Evaluations-es
     ys = [[f(*x) for x in p] for f in fs]