diff --git a/.gitignore b/.gitignore index 570a13e..2aeb440 100644 --- a/.gitignore +++ b/.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 diff --git a/examples/data/README.md b/examples/data/README.md new file mode 100644 index 0000000..4b282d8 --- /dev/null +++ b/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 --crust 1 --dry-run` will tell you how many +blocks it would place before you place any of them. diff --git a/examples/data/moon_apollo15_hires.npz b/examples/data/moon_apollo15_hires.npz deleted file mode 100644 index abb1afd..0000000 Binary files a/examples/data/moon_apollo15_hires.npz and /dev/null differ diff --git a/examples/data/moon_hadley_full.npz b/examples/data/moon_hadley_full.npz deleted file mode 100644 index 4c7d55e..0000000 Binary files a/examples/data/moon_hadley_full.npz and /dev/null differ diff --git a/examples/data/moon_hadley_rille_2m.npz b/examples/data/moon_hadley_rille_2m.npz deleted file mode 100644 index 1dd6a30..0000000 Binary files a/examples/data/moon_hadley_rille_2m.npz and /dev/null differ diff --git a/examples/data/moon_malapert_massif_hires.npz b/examples/data/moon_malapert_massif_hires.npz deleted file mode 100644 index 9009e5e..0000000 Binary files a/examples/data/moon_malapert_massif_hires.npz and /dev/null differ diff --git a/examples/fetch_planet.py b/examples/fetch_planet.py index 3d3ea75..00cd74a 100644 --- a/examples/fetch_planet.py +++ b/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 diff --git a/tests/test_site_dems.py b/tests/test_site_dems.py index 2540384..35f013d 100644 --- a/tests/test_site_dems.py +++ b/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():