# HG changeset patch # User David Barts # Date 1629852000 25200 # Node ID c20b5314774f9152ce5f9d57c0d30bfb74569d2d Initial commit. diff -r 000000000000 -r c20b5314774f Readme.rst --- /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 `_ 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 `_ (specifically, Python 3.9 or better), the +`Pillow Python library `_, and +`GDAL `_ 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. + diff -r 000000000000 -r c20b5314774f make-tile --- /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)