--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Readme.rst	Tue Aug 24 17:40:00 2021 -0700
@@ -0,0 +1,89 @@
+TopoTiler Readme
+================
+
+What Does This Do?
+------------------
+
+This program allows you to turn USGS GeoPDF files (such as those from the
+`National Map <https://nationalmap.gov>`_ program) into tileable GeoTIFF
+files.
+
+By "tileable" I mean suitable for tiling into a complete raster
+layer with gdalwarp or some other such tool. This implies that the so-called
+"map collar" (i.e. the margin and any notation therein, basically anything
+other than the represenataional body of the map) is removed and cropped
+away.
+
+Why Do This?
+------------
+
+It allows one to use such maps as the basis for a raster layer in QGIS
+and other GIS programs.
+
+The National Map products are very nice, high-quality topo maps. Their greatest
+flaw is that they are distributed only as GeoPDF documents, and such a format
+is really only useful for generating printed output (and not as input to a GIS).
+Rendering such documents into tileable raster data is the easiest workaround
+to this drawback.
+
+Why Is This Only a Command-Line Program?
+----------------------------------------
+
+Two reasons:
+
+1. I have not yet had time to write a GUI front end for it. I wrote TopoTiler
+   to serve an occasional need of mine, and a command-line utility works well
+   enough.
+2. I am not yet certain if there is really a need to. A lot of that depends
+   on the need for TopoTiler in the first place. If this program proves useful
+   to a lot of people, and the consensus among its user base is that a GUI
+   would be beneficial, I will probably add one.
+   
+Prerequisites
+-------------
+
+This program requires `python3 <https://python.org/>`_ (specifically, Python 3.9 or better), the
+`Pillow Python library <https://python-pillow.org/>`_, and
+`GDAL <https://gdal.org/>`_ to be installed.
+
+make-tile
+---------
+
+This is the executable script to run.
+
+Synopsis
+^^^^^^^^
+
+``make-tile`` [``-l``/``--layers`` *list*] [``-r``/``--resolution`` *dpi*] *file*
+
+Convert USGS GeoPDF to tileable GeoTIFF. Input files must have an
+extension of ``.pdf`` (case insensitive). Output files will be generated to
+match each input file, with an extension of ``.tiff``.
+
+Arguments
+^^^^^^^^^
+
+-l/--layers
+        Comma-separated list of layers to include. See below for default.
+-r/--resolution
+        Output resolution in DPI. Default is 300.
+        
+Environment
+^^^^^^^^^^^
+
+``GDALBIN`` may be set to the directory containing the GDAL executables. If
+so, they will only be looked for there. Else the standard execution path
+will be searched.
+
+``TMPDIR`` may be set to a place to write temporary files. Doing so is a
+good idea if ``make-tile`` runs out of space in the default temporary files
+area (this program makes very large scratch files).
+
+Layers
+^^^^^^
+
+By default, all ``Map_Frame`` layers, except for
+``Map_Frame.Projection_and_Grids`` and
+``Map_Frame.Terrain.Shaded_Relief``, will be included. This behavior can be
+changed with the ``--layers`` option. 
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/make-tile	Tue Aug 24 17:40:00 2021 -0700
@@ -0,0 +1,266 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# should honor GDALBIN, if set, else use PATH
+# https://opengislab.com/blog/2016/4/2/usgs-geopdf-to-geotif-with-gdal
+# gdalinfo -json -mdd LAYERS gunk.pdf
+# build list of all non-Map.Frame.* layers
+# to that add Map_Frame.Projection_and_Grids
+# First:
+# gdal_translate gunk.pdf /tmp/xxxxx.tif -of GTiff --config GDAL_PDF_LAYERS_OFF Map_Collar,Map_Frame.Projections_and_Grids,Map_Frame.Terrain.Shaded_Relief,Images.Orthoimage --config GDAL_PDF_DPI 300
+# Then:
+# gdal_translate /tmp/xxxxx.tif gunk.tif -of GTiff -srcwin xoff yoff xsize ysize
+# Need to filter out all errors containing: insufficient arguments for Marked Content
+# also fiter out: More than 1000 errors or warnings have been reported. No more will be reported from now.
+
+# I m p o r t s
+
+import os, sys
+import argparse
+import json
+import subprocess, shutil
+import locale
+from abc import abstractmethod
+from enum import IntEnum
+from tempfile import TemporaryDirectory
+
+# PIL is ...weird... and needs both of these.
+import PIL
+from PIL import Image
+
+# V a r i a b l e s
+
+MYNAME = os.path.basename(sys.argv[0])
+estat = 0
+ENCODING = "UTF-8"
+ENVIRON = os.environ.copy()
+for k in ENVIRON.keys():
+    if k.startswith("LC_"): del(ENVIRON[k])
+ENVIRON["LANG"] = ENVIRON["LC_ALL"] = "C." + ENCODING
+
+# C l a s s e s
+
+# A horizontal or vertical strip of an image. Assumes RGB images with 8-bit
+# color.
+class Strip:
+    def __init__(self, based_on, index):
+        if not isinstance(based_on, PIL.Image.Image):
+            raise TypeError("Expecting RGB PIL.Image object for based_on.")
+        if based_on.mode != "RGB":
+            raise ValueError("Expecting RGB PIL.Image object for based_on.")
+        self._index = index
+        self._based_on = based_on
+
+    @property
+    def index(self):
+        return self._index
+
+    @abstractmethod
+    def __len__(self):
+        pass
+
+    @abstractmethod
+    def __getitem__(self, key):
+        pass
+
+    @abstractmethod
+    def __setitem__(self, key, value):
+        pass
+
+    def is_all_white(self):
+        for i in range(len(self)):
+            if self[i] != (255, 255, 255):
+                return False
+        return True
+
+# A horizontal strip of an image.
+class Row(Strip):
+    def __len__(self):
+        return self._based_on.width
+
+    def __getitem__(self, key):
+        return self._based_on.getpixel((key, self._index))
+
+    def __setitem__(self, key, value):
+        self._based_on.setpixel((key, self._index), value)
+
+# A vertical strip of an image.
+class Column(Strip):
+    def __len__(self):
+        return self._based_on.height
+
+    def __getitem__(self, key):
+        return self._based_on.getpixel((self._index, key))
+
+    def __setitem__(self, key, value):
+        self._based_on.setpixel((self._index, key), value)
+
+# A direction in which to search.
+class Direction(IntEnum):
+    DESCENDING = -1
+    ASCENDING = 1
+
+# An orientation in which to search.
+class Orientation(IntEnum):
+    VERTICAL = 0
+    HORIZONTAL = 1
+
+class TopoTilerException(Exception):
+    pass
+
+# F u n c t i o n s
+
+def find_edge(image, orientation, direction):
+    """
+    Find an edge in an image with a white border.
+    """
+    if orientation == Orientation.VERTICAL:
+        what = "column"
+        strip = Column
+        maxn = image.width - 1
+    else:
+        what = "row"
+        strip = Row
+        maxn = image.height - 1
+    # Establish outer bound and sanity-check it.
+    outer = 0 if direction == Direction.ASCENDING else maxn
+    if not strip(image, outer).is_all_white():
+        raise TopoTilerException("Unexpected nonwhite edge ({0} = {1})!".format(what, outer))
+    # Establish inner bound and sanity-check it.
+    inner = outer + direction * maxn // 4
+    if strip(image, inner).is_all_white():
+        raise TopoTilerException("Unexpected white strip ({0} = {1})!".format(what, inner))
+    while abs(outer - inner) > 1:
+        new = (outer + inner) // 2
+        if strip(image, new).is_all_white():
+            outer = new
+        else:
+            inner = new
+    return inner
+
+def positive_int(string):
+    """
+    Specify a positive integer argument.
+    """
+    try:
+        ret = int(string)
+    except ValueError as e:
+        raise argparse.ArgumentTypeError(e.message)
+    if ret <= 0:
+        raise argparse.ArgumentTypeError("%r is 0 or negative" % string)
+    return ret
+
+def gdalcmd(command, *args):
+    """
+    Build gdal command vector for subprocess.Popen
+    """
+    # print('$', command, ' '.join(args))  # debug
+    gdalbin = os.environ.get("GDALBIN")
+    cmd = shutil.which(command) if gdalbin is None else os.path.join(gdalbin, command)
+    ret = list(args)
+    ret.insert(0, cmd)
+    return ret
+
+_MUZZLE = [
+    "insufficient arguments for Marked Content",
+    "More than 1000 errors or warnings have been reported."
+]
+def drainout(stream):
+    """
+    Drain an output/error stream, muzzling the meaningless babble.
+    """
+    for line in stream:
+        printit = True
+        for m in _MUZZLE:
+            if m in line:
+                printit = False
+                break
+        if printit:
+            sys.stderr.write(line)
+
+def waitfor(proc):
+    global estat
+    command_name = os.path.basename(proc.args[0])
+    status = proc.wait()
+    if status < 0:
+        sys.stderr.write("{0}: {1} killed by signal {2}\n".format(MYNAME, command_name, -status))
+        estat = 1
+    elif status > 0:
+        sys.stderr.write("{0}: {1} exited with status {2}\n".format(MYNAME, command_name, status))
+        estat = 1
+
+_BLACKLIST = set(["Map_Frame.Projection_and_Grids", "Map_Frame.Terrain.Shaded_Relief"])
+def get_default_layers(name):
+    """
+    Given a file, analyze it, and get the list of layers to not include
+    when rendering it.
+    """
+    proc = subprocess.Popen(
+        gdalcmd("gdalinfo", "-json", "-mdd", "LAYERS", name),
+        stdout=subprocess.PIPE, stderr=sys.stderr, env=ENVIRON, encoding=ENCODING)
+    data = json.load(proc.stdout)
+    waitfor(proc)
+    layers = set(["Map_Frame"])
+    for layer in data["metadata"]["LAYERS"].values():
+        if layer.startswith("Map_Frame.") and layer not in _BLACKLIST:
+            layers.add(layer)
+    return ",".join(layers)
+
+# M a i n   P r o g r a m
+
+# Parse command-line arguments.
+parser = argparse.ArgumentParser(description='Render GeoPDF into GeoTIFF suitable for tiling.')
+parser.add_argument('--layers', '-l', type=str,
+    help='List of layers to include (see documentation for default).')
+parser.add_argument('--resolution', '-r', type=positive_int, default=300,
+    help='Resolution to render at in DPI.')
+parser.add_argument('file', nargs=1,
+    help='File to read (must be a GeoPDF).')
+try:
+    args = parser.parse_args()
+except SystemExit:
+    sys.exit(2)
+
+# File must end with '.pdf' (case insensitive), or else!
+if not isinstance(args.file, str):
+    # How silly, it gave us a list or a tuple, despite only one of them!
+    args.file = args.file[0]
+if not args.file.lower().endswith('.pdf'):
+    sys.steder.write("{0}: input file must end with .pdf (case insensitive)\n".format(MYNAME))
+
+# Default the set of layers to delete, if necessary
+if args.layers is None:
+    args.layers = get_default_layers(args.file)
+
+with TemporaryDirectory() as td:
+    # Get scratch file name. This goes under TMPDIR; if the default temporary
+    # area is too small, set that environment variable accordingly!
+    tf = os.path.join(td, "temp.tiff")
+
+    # Render. The result will have undesired margins.
+    sys.stdout.write("Rendering (may take a while)...\n")
+    sys.stdout.flush()
+    proc = subprocess.Popen(
+        gdalcmd("gdal_translate", "-q", args.file, tf, "-of", "GTiff", "--config", "GDAL_PDF_LAYERS", args.layers, "--config", "GDAL_PDF_DPI", str(args.resolution)),
+        stdout=subprocess.DEVNULL, stderr=subprocess.PIPE, env=ENVIRON, encoding=ENCODING)
+    drainout(proc.stderr)
+    waitfor(proc)
+    sys.stdout.write("Done!\n")
+
+    # Determine crop marquee.
+    with Image.open(tf) as im:
+        minx = find_edge(im, Orientation.VERTICAL, Direction.ASCENDING)
+        maxx = find_edge(im, Orientation.VERTICAL, Direction.DESCENDING)
+        miny = find_edge(im, Orientation.HORIZONTAL, Direction.ASCENDING)
+        maxy = find_edge(im, Orientation.HORIZONTAL, Direction.DESCENDING)
+
+    # Determine output file name.
+    outf = os.path.splitext(args.file)[0] + ".tiff"
+
+    # Crop.
+    proc = subprocess.Popen(
+        gdalcmd("gdal_translate", "-q", tf, outf, "-of", "GTiff", "-srcwin", str(minx), str(miny), str(maxx-minx+1), str(maxy-miny+1)),
+        stdout=subprocess.DEVNULL, stderr=sys.stderr)
+    waitfor(proc)
+
+# AMF...
+sys.exit(estat)