Follow the strip offsets, and add the other three Apollo DTMs

.gitignore

@@ -32,3 +32,18 @@ examples/data/exo*.csv
 # clear them. Specific to one world and one set of coordinates; useless in a
 # clone and actively wrong if it travelled to another server.
 examples/data/exo_build_*.json
+
+# Landscape windows fetched by fetch_planet.py.
+#
+# The metre-scale ones run to megabytes each and there is no end to how many
+# there could be -- every landing site, every crater, at every window and
+# stride somebody fancies. They are also perfectly reproducible: the source is
+# a fixed published DEM and the reader is deterministic, so the command is the
+# whole definition of the file. examples/data/README.md records the command
+# for each one, which is smaller than the file and easier to read.
+#
+# NOTE, as with the STLs above: this does not untrack anything already
+# committed. The small 118 m and 463 m windows stay in the repository on
+# purpose, so the examples run with no network at all.
+examples/data/moon_*.npz
+examples/data/mars_*.npz

examples/data/README.md

@@ -0,0 +1,108 @@
+# The landscape windows, and how to get them back
+
+`fetch_planet.py` writes an `.npz` here and `build_town.py` reads it. The small
+ones are in the repository; the metre-scale ones are not, because they are
+large, there is no end to how many of them there could be, and git keeps every
+version of everything for ever.
+
+So this file is the record. Every window below is one command, and the command
+is the whole definition of it — same command, same bytes, because the source
+data is a fixed published product and the reader is deterministic.
+
+## Committed
+
+These are small and stay in the repository, so the examples run with no network
+at all.
+
+| file | what | source |
+|---|---|---|
+| `roberts.npz` | Roberts Elementary, Medford MA | USGS lidar + OpenStreetMap |
+| `mars_gale.npz` | Gale Crater | MOLA, 463 m |
+| `mars_olympus.npz` | Olympus Mons | MOLA, 463 m |
+| `mars_hellas.npz` | Hellas Planitia | MOLA, 463 m |
+| `moon_tranquillity.npz` | Tranquility Base | LOLA, 118 m |
+| `moon_apollo15.npz` | Hadley Rille | LOLA, 118 m |
+| `moon_malapert_massif.npz` | Malapert Massif | LOLA, 118 m |
+
+## Fetched on demand
+
+Each is a few seconds and a few megabytes of HTTP range requests. Run these
+from the `examples` directory.
+
+### Apollo 15, Hadley Rille — LROC NAC stereo, 2 m per pixel
+
+```sh
+# the landing site itself: flat, because Apollo 15 landed on a mare plain
+python fetch_planet.py --dem apollo15 --size 300 --name moon_apollo15_hires
+
+# the rille: 323 m deep, and no gaps at this particular window
+python fetch_planet.py --dem apollo15 --lat 26.132 --lon 3.5467 \
+    --size 300 --name moon_hadley_rille_2m
+
+# both at once, coarsened to 4 m so the pair fits in one contiguous map
+python fetch_planet.py --dem apollo15 --lat 26.132 --lon 3.5878 \
+    --size 1480 --rows 740 --stride 2 --name moon_hadley_full
+```
+
+The starting column matters more than it looks. West of about column 400 the
+NAC stereo pair runs out and the coverage becomes ragged, so a window that
+reaches for the far rim of the rille buys it with holes:
+
+| starts at | nodata | width at 4 m |
+|---|---|---|
+| col 420 | 0.02% | 740 blocks |
+| col 380 | 1.45% | 760 |
+| col 340 | 3.97% | 780 |
+| col 300 | 6.37% | 800 |
+
+### Apollo 12, Surveyor Crater — 2 m per pixel
+
+```sh
+python fetch_planet.py --dem apollo12 --lat -3.0157 --lon 336.5847 \
+    --size 600 --name moon_surveyor_crater
+```
+
+208 m across and 21 m rim to floor, measured off the data; the catalogued
+diameter is about 200 m. Conrad put Intrepid on the rim, close enough to walk
+down to Surveyor 3 on the inner slope. None of the six missions landed *in* a
+crater — a crater floor is the last place you want to set down — and this is
+the closest any of them came.
+
+### Apollo 14, Cone Crater — 2 m per pixel
+
+```sh
+python fetch_planet.py --dem apollo14 --lat -3.634361 --lon 342.544006 \
+    --size 1100 --rows 800 --name moon_cone_crater
+```
+
+2.20 x 1.60 km holding both Antares and Cone Crater, 1,481 m apart. The crater
+is 351 m across against a published figure of about 340.
+
+Its depth depends entirely on which rim you stand on, which is worth knowing
+before quoting a number: Cone sits on the flank of a ridge, so it is 44 m from
+the west rim down to the floor, while the ground east of it simply keeps
+climbing and never really has a rim. Measuring from over there gives 88 m,
+which is the ridge and not the crater.
+
+### Malapert Massif — LOLA 5 m south pole site grid
+
+```sh
+python fetch_planet.py --dem malapert --size 300 --name moon_malapert_massif_hires
+```
+
+## Making your own
+
+`--dem` lists what is available; `fetch_planet.py --list` shows the sites and
+the nine further LOLA 5 m grids that are not wired up yet. `--check` prints
+what a DEM covers without reading any pixels, which is the cheap way to find
+out whether the thing you want is in it.
+
+Two switches decide what a window costs:
+
+- `--size` and `--rows` are in PIXELS OF THE SOURCE, not blocks.
+- `--stride` averages n x n of those into one block. A 2 m DEM at `--stride 3`
+  is 6 m to the block, which is how a window covers ground it could not afford
+  at full resolution.
+
+Then `build_town.py --name <name> --crust 1 --dry-run` will tell you how many
+blocks it would place before you place any of them.

examples/fetch_planet.py

@@ -98,6 +98,45 @@
 # GeoTIFF keys say, and unproject() below turns the corner pixels back into
 # latitude and longitude to check they reproduce the published extent.
 SITE_DEMS = {
+    "apollo11": {
+        "body": "moon",
+        "url": "https://pds.lroc.im-ldi.com/data/LRO-L-LROC-5-RDR-V1.0/"
+               "LROLRC_2001/DATA/SDP/NAC_DTM/APOLLO11/NAC_DTM_APOLLO11.TIF",
+        "metres_per_pixel": 2.0,
+        "units_per_metre": 1.0,
+        "lat": 0.674, "lon": 23.473,
+        "what": "Apollo 11, Tranquility Base -- LROC NAC stereo DTM",
+        "source": "LROC / NAC_DTM_APOLLO11, 2 m per pixel",
+        "extent": None,
+        "projection": {"kind": "equirectangular", "radius_m": 1737400.0,
+                       "centre_lon": 180.0, "standard_parallel": 1.0},
+    },
+    "apollo12": {
+        "body": "moon",
+        "url": "https://pds.lroc.im-ldi.com/data/LRO-L-LROC-5-RDR-V1.0/"
+               "LROLRC_2001/DATA/SDP/NAC_DTM/APOLLO12/NAC_DTM_APOLLO12.TIF",
+        "metres_per_pixel": 2.0,
+        "units_per_metre": 1.0,
+        "lat": -3.012, "lon": 336.578,
+        "what": "Apollo 12, Surveyor Crater -- LROC NAC stereo DTM",
+        "source": "LROC / NAC_DTM_APOLLO12, 2 m per pixel",
+        "extent": None,
+        "projection": {"kind": "equirectangular", "radius_m": 1737400.0,
+                       "centre_lon": 180.0, "standard_parallel": -3.0},
+    },
+    "apollo14": {
+        "body": "moon",
+        "url": "https://pds.lroc.im-ldi.com/data/LRO-L-LROC-5-RDR-V1.0/"
+               "LROLRC_2001/DATA/SDP/NAC_DTM/APOLLO14/NAC_DTM_APOLLO14.TIF",
+        "metres_per_pixel": 2.0,
+        "units_per_metre": 1.0,
+        "lat": -3.645, "lon": 342.522,
+        "what": "Apollo 14, Fra Mauro and Cone Crater -- LROC NAC stereo DTM",
+        "source": "LROC / NAC_DTM_APOLLO14, 2 m per pixel",
+        "extent": None,
+        "projection": {"kind": "equirectangular", "radius_m": 1737400.0,
+                       "centre_lon": 180.0, "standard_parallel": -3.0},
+    },
     "apollo15": {
         "body": "moon",
         "url": "https://pds.lroc.im-ldi.com/data/LRO-L-LROC-5-RDR-V1.0/"
@@ -318,14 +357,46 @@ def read_window(url, info, row0, col0, rows, cols, wrap=True):
     """
     bo = info["byte_order"]
     width = info["width"]
-    row_bytes = info["strip_bytes"][0]
-
-    start = info["strip_offsets"][row0]
-    end = info["strip_offsets"][row0 + rows - 1] + row_bytes - 1
-    raw = _get(url, start, end)
-
+    offsets = info["strip_offsets"]
+    counts = info["strip_bytes"]
     dtype = np.dtype(bo + info["sample_dtype"])
-    block = np.frombuffer(raw, dtype=dtype, count=rows * width).reshape(rows, width)
+
+    # Rows are NOT necessarily in file order.
+    #
+    # Apollo 15 stores them one after another, which is what this used to
+    # assume: take the offset of the first row, the offset of the last, and
+    # ask for everything between. Apollo 12 keeps row 0 at the very END of the
+    # file, after every other row, with a gap in the middle for good measure.
+    # Reading it as one run returned bytes that were all real numbers from
+    # somewhere in the image, just not the rows asked for -- elevations came
+    # out at 3e38 and the landscape was noise.
+    #
+    # So: group the wanted rows into runs that really are contiguous, and ask
+    # for one byte range per run. A file that is in order still costs exactly
+    # one request. Runs are also capped, because a single very large range
+    # comes back truncated from this server and numpy then reads off the end.
+    MAX_RUN = 8 << 20
+    runs = []
+    for r in range(row0, row0 + rows):
+        off, cnt = offsets[r], counts[r]
+        if runs and off == runs[-1][1] and (off + cnt - runs[-1][0]) <= MAX_RUN:
+            runs[-1][1] = off + cnt
+            runs[-1][2].append(r)
+        else:
+            runs.append([off, off + cnt, [r]])
+
+    rowdata = {}
+    for start, end, wanted in runs:
+        raw = _get(url, start, end - 1)
+        if len(raw) < end - start:
+            raise ValueError("short read: asked for %d bytes, got %d"
+                             % (end - start, len(raw)))
+        pos = 0
+        for r in wanted:
+            rowdata[r] = np.frombuffer(raw, dtype=dtype, count=width, offset=pos)
+            pos += counts[r]
+
+    block = np.stack([rowdata[r] for r in range(row0, row0 + rows)])
     if wrap:
         # A global mosaic joins up, so a window across the antimeridian works.
         idx = np.arange(col0, col0 + cols) % width

tests/test_site_dems.py

@@ -42,6 +42,23 @@
 }
 MALAPERT_PROJ = fetch_planet.SITE_DEMS["malapert"]["projection"]
 
+# The other three Apollo DTMs, again straight from the files. Apollo 12's is
+# the odd one: a compact 4.94 x 2.73 km tile rather than a long strip, and the
+# only one whose rows are not stored in order.
+HEADERS = {
+    "apollo15": APOLLO15,
+    "malapert": MALAPERT,
+    "apollo11": {"width": 2111, "height": 13978,
+                 "x_min": -4748754.0000014, "y_max": 37496.000000011,
+                 "x_scale": 2.0000000000006, "y_scale": 2.0000000000006},
+    "apollo12": {"width": 2468, "height": 1367,
+                 "x_min": 4737648.0000014, "y_max": -89822.000000027,
+                 "x_scale": 2.0000000000006, "y_scale": 2.0000000000006},
+    "apollo14": {"width": 2156, "height": 14400,
+                 "x_min": 4919932.0000015, "y_max": -89856.000000027,
+                 "x_scale": 2.0000000000006, "y_scale": 2.0000000000006},
+}
+
 
 def _full(info):
     """Fill in the derived extent the reader normally computes."""
@@ -95,12 +112,13 @@ def test_projection_round_trips(name, info, proj, lat, lon):
 
 def test_the_landing_site_is_in_its_own_dem():
     """Each site DEM must actually contain the coordinates it is filed under."""
+    assert set(HEADERS) == set(fetch_planet.SITE_DEMS),         "a DEM was added without its header, so nothing here checks it"
     for name, dem in fetch_planet.SITE_DEMS.items():
-        info = _full(APOLLO15 if name == "apollo15" else MALAPERT)
+        info = _full(HEADERS[name])
         col, row = fetch_planet.pixel_for(dem["lat"], dem["lon"], info,
                                           dem["projection"])
-        assert 0 <= col < info["width"], name
-        assert 0 <= row < info["height"], name
+        assert 0 <= col < info["width"], "%s: col %d of %d" % (name, col, info["width"])
+        assert 0 <= row < info["height"], "%s: row %d of %d" % (name, row, info["height"])
 
 
 def test_polar_stereographic_pole_is_the_origin():