view make-tile @ 8:290b0b2a584b

Fix bug.
author David Barts <n5jrn@me.com>
date Wed, 25 Aug 2021 08:33:59 -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)