Mercurial > cgi-bin > hgweb.cgi > TopoTiler
view make-tile @ 3:ab4d4434c37f
This is useful when grabbing coordinates from QGIS.
author | David Barts <n5jrn@me.com> |
---|---|
date | Wed, 25 Aug 2021 00:24:30 -0700 |
parents | 94762476b171 |
children | 1944acce0e6f |
line wrap: on
line source
#!/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): """ Wait for a GDAL command to finish. """ 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 (default: 300).') 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)