Coarser windows, a crust that stays up, and getBlocks in the right order

examples/build_town.py

@@ -61,6 +61,43 @@
     },
 }
 
+# Blocks that fall when nothing holds them up, and what to use instead.
+#
+# This matters only because a crust floats, and it very nearly does not
+# announce itself. Blocks are placed with physics off -- deliberately, because
+# running block physics for every one of half a million placements is what
+# stalled the server and got it killed by the watchdog -- so gravity is never
+# evaluated at build time and a crust of regolith sits in mid-air looking
+# perfectly solid. It stays that way until anything at all causes a block
+# update nearby. Mine one block and the update spreads, every neighbour
+# discovers it is unsupported, and the landscape drains away.
+#
+# The Moon's top three layers were concrete powder and Mars's top two were red
+# sand, so both were a landslide waiting for the first pickaxe. The substitutes
+# are chosen to keep the colour: the point of the palette is that regolith is
+# grey, not that it is dust.
+_COLOURS = ("WHITE", "ORANGE", "MAGENTA", "LIGHT_BLUE", "YELLOW", "LIME",
+            "PINK", "GRAY", "LIGHT_GRAY", "CYAN", "PURPLE", "BLUE", "BROWN",
+            "GREEN", "RED", "BLACK")
+FALLING = {"%s_CONCRETE_POWDER" % c: "%s_CONCRETE" % c for c in _COLOURS}
+FALLING.update({
+    "SAND": "SANDSTONE",
+    "SUSPICIOUS_SAND": "SANDSTONE",
+    "RED_SAND": "RED_SANDSTONE",
+    "GRAVEL": "COBBLESTONE",
+    "SUSPICIOUS_GRAVEL": "COBBLESTONE",
+    "ANVIL": "IRON_BLOCK",
+    "POINTED_DRIPSTONE": "DRIPSTONE_BLOCK",
+    "SCAFFOLDING": "OAK_PLANKS",
+    "DRAGON_EGG": "OBSIDIAN",
+})
+
+
+def anchored(material):
+    """The nearest block of the same colour that stays where it is put."""
+    return FALLING.get(material, material)
+
+
 WALL = "WHITE_CONCRETE"
 ROOF = "GRAY_CONCRETE"
 
@@ -347,17 +384,51 @@ def to_voxels(ground, buildings, base_y, depth, hollow=True, body="earth",
     # surface either way; the difference is only visible from underneath, and
     # from underneath there is nothing to see.
     #
+    # But a FLAT crust is both too much and too little. Too much because on
+    # level ground one block is plenty -- nothing can see past it. Too little
+    # because where the ground drops sharply, a column whose neighbour sits
+    # more than `crust` blocks below it has daylight under its lip, and you
+    # see straight through the wall of the valley you came to look at.
+    #
+    # So `crust` is a MINIMUM, and each column is thickened to cover the drop
+    # to its lowest neighbour. On Hadley Rille that turns 1.37M blocks into
+    # 310k -- four times cheaper -- while closing 671 columns of holes that a
+    # flat crust of five left open, because the deepest drop there is 81
+    # blocks and no sane constant covers it.
+    #
     # The label is i + 1, not len(soil) - i. The latter reads the list
     # backwards and buries the grass: yards came out coarse dirt with the turf
     # three blocks down, and Mars had terracotta on top instead of red sand.
     # Anywhere with a land-cover tag was painted over afterwards and looked
     # right, so only the untagged ground -- which is to say people's gardens --
     # showed it.
+    # How far below each column's top its lowest neighbour sits. Off the edge
+    # of the map there is no neighbour, so nothing is required there and the
+    # rim comes out thin, which is what a cut-out of terrain should look like.
+    thickness = None
+    if crust is not None:
+        drop = np.zeros_like(ground_h)
+        for axis in (0, 1):
+            for step in (1, -1):
+                shifted = np.roll(ground_h, step, axis=axis)
+                if axis == 0:
+                    if step == 1:
+                        shifted[0, :] = ground_h[0, :]
+                    else:
+                        shifted[-1, :] = ground_h[-1, :]
+                else:
+                    if step == 1:
+                        shifted[:, 0] = ground_h[:, 0]
+                    else:
+                        shifted[:, -1] = ground_h[:, -1]
+                drop = np.maximum(drop, ground_h - shifted)
+        thickness = np.maximum(int(crust), np.maximum(drop, 0))
+
     for x in range(nx):
         col = ground_h[x]
         for z in range(nz):
             top = col[z]
-            bottom = 0 if crust is None else max(0, top - int(crust))
+            bottom = 0 if crust is None else max(0, top - int(thickness[x, z]))
             world[x, bottom:top, z] = stone
             for i, _ in enumerate(soil):
                 y = top - 1 - i
@@ -537,6 +608,11 @@ def label_for(material):
                             world[px, py, pz] = leaf_label
 
     palette = soil + [deep, WALL, ROOF] + extra
+    if crust is not None:
+        # Only when there is a crust. A build filled to the floor rests on
+        # something, so sand on Mars can behave like sand, which is half the
+        # reason it is sand.
+        palette = [anchored(m) for m in palette]
     return world, palette
 
 

examples/exoplanets.py

@@ -102,7 +102,10 @@
 ROCK_BY_TEMP = (
     (180.0, ["PACKED_ICE", "BLUE_ICE"]),
     (320.0, ["GRASS_BLOCK", "STONE"]),
-    (600.0, ["TERRACOTTA", "RED_SAND"]),
+    # Red sandstone, not red sand: these hang in the sky with nothing under
+    # them, and a gravity-affected block up there is one block update away
+    # from raining the planet onto the ground.
+    (600.0, ["TERRACOTTA", "RED_SANDSTONE"]),
     (1200.0, ["BLACKSTONE", "BASALT"]),
     (float("inf"), ["MAGMA_BLOCK", "NETHERRACK"]),
 )

examples/fetch_planet.py

@@ -45,6 +45,7 @@
 import sys
 import urllib.error
 import urllib.request
+import warnings
 
 try:
     import numpy as np
@@ -281,6 +282,32 @@ def at(off, n):
     return info
 
 
+def downsample(grid, stride):
+    """Average stride x stride pixels into one, ignoring gaps.
+
+    Taking every nth pixel instead would be quicker and wrong: it is point
+    sampling a surface with real detail in it, so a narrow ridge either
+    survives at full height or vanishes entirely depending on where the grid
+    happens to fall. Averaging keeps the shape and quietly loses the detail,
+    which is what a coarser map is supposed to do.
+
+    A window is trimmed to a whole number of blocks first, so the last part
+    row is dropped rather than averaged against nothing.
+    """
+    stride = int(stride)
+    if stride <= 1:
+        return grid
+    rows = (grid.shape[0] // stride) * stride
+    cols = (grid.shape[1] // stride) * stride
+    block = grid[:rows, :cols].reshape(rows // stride, stride,
+                                       cols // stride, stride)
+    with warnings.catch_warnings():
+        # A block that is entirely nodata averages to nan, which is correct
+        # and which numpy would rather warn about every time.
+        warnings.simplefilter("ignore", RuntimeWarning)
+        return np.nanmean(block, axis=(1, 3))
+
+
 def read_window(url, info, row0, col0, rows, cols, wrap=True):
     """A rows x cols window, read as one contiguous range request.
 
@@ -409,7 +436,16 @@ def main():
     p.add_argument("--place", help="a named landing site or feature; --list to see them")
     p.add_argument("--lat", type=float, help="centre latitude, if not using --place")
     p.add_argument("--lon", type=float, help="centre longitude, degrees east")
-    p.add_argument("--size", type=int, default=400, help="window across, in pixels (default 400)")
+    p.add_argument("--size", type=int, default=400,
+                   help="window across, in PIXELS of the source map (default 400)")
+    p.add_argument("--rows", type=int,
+                   help="window north-south, in pixels (default: same as --size). "
+                        "The interesting axis is often only one of them")
+    p.add_argument("--stride", type=int, default=1,
+                   help="average this many pixels square 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 (default 1)")
     p.add_argument("--vscale", type=float, default=1.0,
                    help="vertical exaggeration; 1 is true scale (default 1)")
     p.add_argument("--name", help="output name (default: body_place)")
@@ -491,8 +527,10 @@ def main():
     if args.check:
         return
 
+    rows_wanted = args.rows or args.size
     cx, cy = pixel_for(lat, lon, info, proj)
     half = args.size // 2
+    half_rows = rows_wanted // 2
 
     # A global mosaic wraps and is huge, so a window always fits. A site DEM
     # is a few thousand pixels across and asking for the middle of a 21 km
@@ -503,38 +541,57 @@ def main():
             sys.exit("%.4f, %.4f is not inside %s (pixel %d,%d of %d x %d).\n"
                      "Use --check to see what it covers."
                      % (lat, lon, args.dem, cx, cy, info["width"], info["height"]))
-        if args.size > min(info["width"], info["height"]):
-            print("  note: --size %d is larger than the DEM (%d x %d); trimming"
-                  % (args.size, info["width"], info["height"]))
-            args.size = min(info["width"], info["height"])
+        if args.size > info["width"]:
+            print("  note: --size %d is wider than the DEM (%d); trimming"
+                  % (args.size, info["width"]))
+            args.size = info["width"]
             half = args.size // 2
+        if rows_wanted > info["height"]:
+            rows_wanted = info["height"]
+            half_rows = rows_wanted // 2
         col0 = max(0, min(cx - half, info["width"] - args.size))
     else:
         col0 = cx - half
-    row0 = max(0, min(cy - half, info["height"] - args.size))
-
-    span_km = args.size * mpp / 1000.0
-    mb = args.size * info["strip_bytes"][0] / 1e6
-    print("  window %d x %d pixels = %.0f x %.0f km  (about %.0f MB of range requests)"
-          % (args.size, args.size, span_km, span_km, mb))
-
-    raw = read_window(url, info, row0, col0, args.size, args.size, wrap=dem is None)
-    ground = raw / source["units_per_metre"]
-
-    # Gaps. The global mosaics use the most negative 16-bit integer; the site
-    # DEMs say what they use in the TIFF, and one of them says "nan", which
-    # never equals itself and so has to be tested for separately.
+    row0 = max(0, min(cy - half_rows, info["height"] - rows_wanted))
+
+    span_x = args.size * mpp / 1000.0
+    span_z = rows_wanted * mpp / 1000.0
+    mb = rows_wanted * info["strip_bytes"][0] / 1e6
+    print("  window %d x %d pixels = %.2f x %.2f km  (about %.0f MB of range requests)"
+          % (args.size, rows_wanted, span_x, span_z, mb))
+
+    raw = read_window(url, info, row0, col0, rows_wanted, args.size,
+                      wrap=dem is None)
+
+    # Gaps become NaN BEFORE anything averages them. The global mosaics mark
+    # nodata with the most negative 16-bit integer; the site DEMs say what
+    # they use in the TIFF, and for the NAC DTM that is -3.4e38.
+    #
+    # Averaging first is a quiet disaster: one nodata pixel in a 2x2 block
+    # drags the average to -8.5e37, which is not obviously a gap, is not
+    # caught by a later "is this the nodata value" test because it no longer
+    # equals it, and builds as a hole a hundred billion times deeper than the
+    # crater. NaN spreads harmlessly instead, and nanmean ignores it.
     nodata = info.get("nodata")
     if nodata is not None and nodata == nodata:
-        ground[raw <= nodata + abs(nodata) * 1e-6] = np.nan
-    ground[~np.isfinite(raw)] = np.nan
+        raw[raw <= nodata + abs(nodata) * 1e-6] = np.nan
+    raw[~np.isfinite(raw)] = np.nan
+
+    if args.stride > 1:
+        before = raw.shape
+        raw = downsample(raw, args.stride)
+        mpp = mpp * args.stride
+        print("  averaged %dx%d pixels per block: %s -> %s at %g m per block"
+              % (args.stride, args.stride, "x".join(map(str, before)),
+                 "x".join(map(str, raw.shape)), mpp))
+    ground = raw / source["units_per_metre"]
     gaps = int(np.isnan(ground).sum())
     if gaps:
         print("  %d of %d pixels have no data (%.1f%%)"
               % (gaps, ground.size, 100.0 * gaps / ground.size))
     lo, hi = float(np.nanmin(ground)), float(np.nanmax(ground))
-    print("  elevation %.0f to %.0f m  (relief %.0f m over %.0f km)"
-          % (lo, hi, hi - lo, span_km))
+    print("  elevation %.0f to %.0f m  (relief %.0f m over %.2f km)"
+          % (lo, hi, hi - lo, max(span_x, span_z)))
 
     if args.vscale != 1.0:
         mid = (lo + hi) / 2.0
@@ -556,6 +613,7 @@ def main():
         meta=json.dumps({
             "body": args.body, "place": args.place, "what": what,
             "lat": lat, "lon": lon, "size": args.size,
+            "rows": rows_wanted, "stride": args.stride,
             "metres_per_pixel": mpp, "vscale": args.vscale,
             "elevation_min_m": lo, "elevation_max_m": hi,
             "metres_per_cell": mpp,

examples/moon_landings.py

@@ -111,12 +111,9 @@ def build_apollo_lm(scale=0.5):
 
     # Four landing legs, out and down, with a round footpad on the end.
     for sx, sz in ((1, 1), (1, -1), (-1, 1), (-1, -1)):
-        steps = max(leg_span - desc_r, 1)
-        for i in range(steps + 1):
-            t = i / steps
-            x = int(round(sx * (desc_r + (leg_span - desc_r) * t)))
-            z = int(round(sz * (desc_r + (leg_span - desc_r) * t)))
-            y = int(round(desc_h * 0.55 * (1 - t)))
+        top = (sx * desc_r, int(round(desc_h * 0.55)), sz * desc_r)
+        foot = (sx * leg_span, 0, sz * leg_span)
+        for x, y, z in connected(top, foot):
             out.append((x, y, z, "IRON_BARS"))
         px, pz = sx * leg_span, sz * leg_span
         for dx in (-1, 0, 1):
@@ -152,6 +149,27 @@ def build_apollo_lm(scale=0.5):
     return out
 
 
+def connected(a, b):
+    """Every point from a to b, each sharing a FACE with the one before it.
+
+    Sweeping a parametric line and rounding gives points that step diagonally,
+    and two blocks meeting only at an edge are not neighbours as far as
+    Minecraft is concerned. Iron bars work out what to join onto by looking at
+    their four side faces, so a landing leg drawn that way came out as a row
+    of loose rods hanging in the air with daylight between them.
+
+    Moving one axis at a time -- the one with furthest still to go -- costs a
+    few more blocks and makes the leg look like a leg.
+    """
+    cur = list(a)
+    pts = [tuple(cur)]
+    while tuple(cur) != tuple(b):
+        far = max(range(3), key=lambda i: abs(b[i] - cur[i]))
+        cur[far] += 1 if b[far] > cur[far] else -1
+        pts.append(tuple(cur))
+    return pts
+
+
 def build_lrv(scale=0.5):
     """The Lunar Roving Vehicle: two seats, four wire wheels, a dish aerial."""
     b = lambda m: max(int(round(m / scale)), 1)

pyncraft/minecraft.py

@@ -194,6 +194,28 @@ def pollChatPosts(self):
         return [ChatEvent.Post(int(e[:e.find(",")]), e[e.find(",") + 1:]) for e in events]
 
 
+def reshape_blocks(names, xs, ys, zs):
+    """The server's flat block list into [y][x][z], which is the order it uses.
+
+    Kept out of the Minecraft class so it can be tested without a server: the
+    whole bug was in the arithmetic, and the arithmetic needs no socket.
+    """
+    if len(names) != xs * ys * zs:
+        # `warn` here is this package's logger, not warnings.warn: the
+        # `from .logger import *` above shadows it. That is worth knowing
+        # before writing a test that waits for a warning which never arrives.
+        warn("getBlocks: expected %d blocks, got %d" % (xs * ys * zs, len(names)))
+    out, k = [], 0
+    for _ in range(ys):
+        slab = []
+        for _ in range(xs):
+            slab.append(names[k:k + zs])
+            k += zs
+        out.append(slab)
+    return out
+
+
+
 class Minecraft:
     """The main class to interact with a running instance of Minecraft Pi."""
     def __init__(self, connection, playerId):
@@ -259,25 +281,23 @@ def getBlockData(self, x, y=None, z=None, parse: bool = True) -> dict:
         return data
 
     def getBlocks(self, x1:int, y1:int, z1:int, x2:int, y2:int, z2:int) -> list:
-        """Get a cuboid of blocks (x0,y0,z0,x1,y1,z1) => [id:int]"""
+        """A cuboid of blocks, as nested lists indexed [y][x][z].
+
+        The order is the server's, and it is not the obvious one: FruitJuice
+        walks the cuboid with Y outermost, then X, then Z. So the flat reply is
+        a stack of horizontal slabs, each slab a set of rows running east, each
+        row running south.
+
+        This used to reshape into slabs of xSize * ySize, which is the right
+        SIZE only when the cuboid happens to be as tall as it is deep. Any
+        other shape silently came back scrambled -- every block was a real
+        block from somewhere in the box, just not from where you thought, so
+        nothing looked wrong until the answer was checked against a place
+        somebody could actually stand.
+        """
         blocks = self.conn.sendReceive(b"world.getBlocks", x1, y1, z1, x2, y2, z2)
-        arr1d = blocks.split(',')
-
-        xSize = abs(x1 - x2) + 1
-        ySize = abs(y1 - y2) + 1
-        zSize = abs(z1 - z2) + 1
-        totalSize = xSize * ySize * zSize
-        arr3d = []
-
-        if len(arr1d) != totalSize:
-            warn('Get number of blocks is incomplete')
-
-        for i in range(0,totalSize,xSize*ySize):
-            curArr = []
-            for j in range(0,xSize*ySize,xSize):
-                curArr.append(arr1d[i+j:i+j+xSize])
-            arr3d.append(curArr)
-        return arr3d
+        return reshape_blocks(blocks.split(","),
+                              abs(x1 - x2) + 1, abs(y1 - y2) + 1, abs(z1 - z2) + 1)
 
     # DIRECTION: NORTH SOUTH EAST WEST
     # FACE: FLOOR, CEILING, WALL