diff --git a/README.md b/README.md index 16d2168..7c3f9a6 100644 --- a/README.md +++ b/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. diff --git a/interpolation.py b/interpolation.py index ae2b678..ff3e735 100755 --- a/interpolation.py +++ b/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,7 +141,7 @@ 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] @@ -157,7 +149,7 @@ 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] @@ -165,8 +157,35 @@ 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') diff --git a/nevis/_bng.py b/nevis/_bng.py index b530d00..f0b2dc4 100644 --- a/nevis/_bng.py +++ b/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), } diff --git a/nevis/_nevis_version.py b/nevis/_nevis_version.py index 8bc951c..524a574 100644 --- a/nevis/_nevis_version.py +++ b/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__) diff --git a/nevis/_os_terrain_50.py b/nevis/_os_terrain_50.py index 2061d65..24ef94c 100644 --- a/nevis/_os_terrain_50.py +++ b/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): diff --git a/nevis/_plot.py b/nevis/_plot.py index 0610ed8..7da687d 100644 --- a/nevis/_plot.py +++ b/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]