|
0
|
1 #!/usr/bin/env python3
|
|
|
2 # -*- coding: utf-8 -*-
|
|
12
|
3 # Make a tileable GeoTIFF from a USGS National Map GeoPDF.
|
|
0
|
4
|
|
|
5 # I m p o r t s
|
|
|
6
|
|
|
7 import os, sys
|
|
|
8 import argparse
|
|
|
9 import json
|
|
|
10 import subprocess, shutil
|
|
|
11 import locale
|
|
|
12 from abc import abstractmethod
|
|
|
13 from enum import IntEnum
|
|
|
14 from tempfile import TemporaryDirectory
|
|
|
15
|
|
|
16 # PIL is ...weird... and needs both of these.
|
|
|
17 import PIL
|
|
|
18 from PIL import Image
|
|
|
19
|
|
|
20 # V a r i a b l e s
|
|
|
21
|
|
|
22 MYNAME = os.path.basename(sys.argv[0])
|
|
|
23 estat = 0
|
|
|
24 ENCODING = "UTF-8"
|
|
|
25 ENVIRON = os.environ.copy()
|
|
|
26 for k in ENVIRON.keys():
|
|
|
27 if k.startswith("LC_"): del(ENVIRON[k])
|
|
|
28 ENVIRON["LANG"] = ENVIRON["LC_ALL"] = "C." + ENCODING
|
|
|
29
|
|
|
30 # C l a s s e s
|
|
|
31
|
|
|
32 # A horizontal or vertical strip of an image. Assumes RGB images with 8-bit
|
|
|
33 # color.
|
|
|
34 class Strip:
|
|
|
35 def __init__(self, based_on, index):
|
|
|
36 if not isinstance(based_on, PIL.Image.Image):
|
|
|
37 raise TypeError("Expecting RGB PIL.Image object for based_on.")
|
|
|
38 if based_on.mode != "RGB":
|
|
|
39 raise ValueError("Expecting RGB PIL.Image object for based_on.")
|
|
|
40 self._index = index
|
|
|
41 self._based_on = based_on
|
|
|
42
|
|
|
43 @property
|
|
|
44 def index(self):
|
|
|
45 return self._index
|
|
|
46
|
|
|
47 @abstractmethod
|
|
|
48 def __len__(self):
|
|
|
49 pass
|
|
|
50
|
|
|
51 @abstractmethod
|
|
|
52 def __getitem__(self, key):
|
|
|
53 pass
|
|
|
54
|
|
|
55 @abstractmethod
|
|
|
56 def __setitem__(self, key, value):
|
|
|
57 pass
|
|
|
58
|
|
|
59 def is_all_white(self):
|
|
|
60 for i in range(len(self)):
|
|
|
61 if self[i] != (255, 255, 255):
|
|
|
62 return False
|
|
|
63 return True
|
|
|
64
|
|
|
65 # A horizontal strip of an image.
|
|
|
66 class Row(Strip):
|
|
|
67 def __len__(self):
|
|
|
68 return self._based_on.width
|
|
|
69
|
|
|
70 def __getitem__(self, key):
|
|
|
71 return self._based_on.getpixel((key, self._index))
|
|
|
72
|
|
|
73 def __setitem__(self, key, value):
|
|
|
74 self._based_on.setpixel((key, self._index), value)
|
|
|
75
|
|
|
76 # A vertical strip of an image.
|
|
|
77 class Column(Strip):
|
|
|
78 def __len__(self):
|
|
|
79 return self._based_on.height
|
|
|
80
|
|
|
81 def __getitem__(self, key):
|
|
|
82 return self._based_on.getpixel((self._index, key))
|
|
|
83
|
|
|
84 def __setitem__(self, key, value):
|
|
|
85 self._based_on.setpixel((self._index, key), value)
|
|
|
86
|
|
|
87 # A direction in which to search.
|
|
|
88 class Direction(IntEnum):
|
|
|
89 DESCENDING = -1
|
|
|
90 ASCENDING = 1
|
|
|
91
|
|
|
92 # An orientation in which to search.
|
|
|
93 class Orientation(IntEnum):
|
|
|
94 VERTICAL = 0
|
|
|
95 HORIZONTAL = 1
|
|
|
96
|
|
|
97 class TopoTilerException(Exception):
|
|
|
98 pass
|
|
|
99
|
|
|
100 # F u n c t i o n s
|
|
|
101
|
|
|
102 def find_edge(image, orientation, direction):
|
|
|
103 """
|
|
|
104 Find an edge in an image with a white border.
|
|
|
105 """
|
|
|
106 if orientation == Orientation.VERTICAL:
|
|
|
107 what = "column"
|
|
|
108 strip = Column
|
|
11
|
109 maxn = image.width
|
|
0
|
110 else:
|
|
|
111 what = "row"
|
|
|
112 strip = Row
|
|
11
|
113 maxn = image.height
|
|
0
|
114 # Establish outer bound and sanity-check it.
|
|
11
|
115 outer = 0 if direction == Direction.ASCENDING else maxn - 1
|
|
0
|
116 if not strip(image, outer).is_all_white():
|
|
|
117 raise TopoTilerException("Unexpected nonwhite edge ({0} = {1})!".format(what, outer))
|
|
|
118 # Establish inner bound and sanity-check it.
|
|
|
119 inner = outer + direction * maxn // 4
|
|
|
120 if strip(image, inner).is_all_white():
|
|
|
121 raise TopoTilerException("Unexpected white strip ({0} = {1})!".format(what, inner))
|
|
|
122 while abs(outer - inner) > 1:
|
|
|
123 new = (outer + inner) // 2
|
|
|
124 if strip(image, new).is_all_white():
|
|
|
125 outer = new
|
|
|
126 else:
|
|
|
127 inner = new
|
|
|
128 return inner
|
|
|
129
|
|
|
130 def positive_int(string):
|
|
|
131 """
|
|
|
132 Specify a positive integer argument.
|
|
|
133 """
|
|
|
134 try:
|
|
|
135 ret = int(string)
|
|
|
136 except ValueError as e:
|
|
|
137 raise argparse.ArgumentTypeError(e.message)
|
|
|
138 if ret <= 0:
|
|
|
139 raise argparse.ArgumentTypeError("%r is 0 or negative" % string)
|
|
|
140 return ret
|
|
|
141
|
|
|
142 def gdalcmd(command, *args):
|
|
|
143 """
|
|
|
144 Build gdal command vector for subprocess.Popen
|
|
|
145 """
|
|
|
146 # print('$', command, ' '.join(args)) # debug
|
|
|
147 gdalbin = os.environ.get("GDALBIN")
|
|
|
148 cmd = shutil.which(command) if gdalbin is None else os.path.join(gdalbin, command)
|
|
|
149 ret = list(args)
|
|
|
150 ret.insert(0, cmd)
|
|
|
151 return ret
|
|
|
152
|
|
|
153 _MUZZLE = [
|
|
|
154 "insufficient arguments for Marked Content",
|
|
|
155 "More than 1000 errors or warnings have been reported."
|
|
|
156 ]
|
|
|
157 def drainout(stream):
|
|
|
158 """
|
|
|
159 Drain an output/error stream, muzzling the meaningless babble.
|
|
|
160 """
|
|
|
161 for line in stream:
|
|
|
162 printit = True
|
|
|
163 for m in _MUZZLE:
|
|
|
164 if m in line:
|
|
|
165 printit = False
|
|
|
166 break
|
|
|
167 if printit:
|
|
|
168 sys.stderr.write(line)
|
|
|
169
|
|
|
170 def waitfor(proc):
|
|
1
|
171 """
|
|
|
172 Wait for a GDAL command to finish.
|
|
|
173 """
|
|
0
|
174 global estat
|
|
|
175 command_name = os.path.basename(proc.args[0])
|
|
|
176 status = proc.wait()
|
|
|
177 if status < 0:
|
|
|
178 sys.stderr.write("{0}: {1} killed by signal {2}\n".format(MYNAME, command_name, -status))
|
|
|
179 estat = 1
|
|
|
180 elif status > 0:
|
|
|
181 sys.stderr.write("{0}: {1} exited with status {2}\n".format(MYNAME, command_name, status))
|
|
|
182 estat = 1
|
|
|
183
|
|
|
184 _BLACKLIST = set(["Map_Frame.Projection_and_Grids", "Map_Frame.Terrain.Shaded_Relief"])
|
|
|
185 def get_default_layers(name):
|
|
|
186 """
|
|
|
187 Given a file, analyze it, and get the list of layers to not include
|
|
|
188 when rendering it.
|
|
|
189 """
|
|
|
190 proc = subprocess.Popen(
|
|
|
191 gdalcmd("gdalinfo", "-json", "-mdd", "LAYERS", name),
|
|
|
192 stdout=subprocess.PIPE, stderr=sys.stderr, env=ENVIRON, encoding=ENCODING)
|
|
|
193 data = json.load(proc.stdout)
|
|
|
194 waitfor(proc)
|
|
|
195 layers = set(["Map_Frame"])
|
|
|
196 for layer in data["metadata"]["LAYERS"].values():
|
|
|
197 if layer.startswith("Map_Frame.") and layer not in _BLACKLIST:
|
|
|
198 layers.add(layer)
|
|
|
199 return ",".join(layers)
|
|
|
200
|
|
|
201 # M a i n P r o g r a m
|
|
|
202
|
|
|
203 # Parse command-line arguments.
|
|
|
204 parser = argparse.ArgumentParser(description='Render GeoPDF into GeoTIFF suitable for tiling.')
|
|
|
205 parser.add_argument('--layers', '-l', type=str,
|
|
|
206 help='List of layers to include (see documentation for default).')
|
|
10
|
207 parser.add_argument('--output', '-o', type=str,
|
|
|
208 help='Name of output file (default: *.tiff).')
|
|
0
|
209 parser.add_argument('--resolution', '-r', type=positive_int, default=300,
|
|
2
|
210 help='Resolution to render at in DPI (default: 300).')
|
|
0
|
211 parser.add_argument('file', nargs=1,
|
|
|
212 help='File to read (must be a GeoPDF).')
|
|
|
213 try:
|
|
|
214 args = parser.parse_args()
|
|
|
215 except SystemExit:
|
|
|
216 sys.exit(2)
|
|
|
217
|
|
|
218 # File must end with '.pdf' (case insensitive), or else!
|
|
|
219 if not isinstance(args.file, str):
|
|
|
220 # How silly, it gave us a list or a tuple, despite only one of them!
|
|
|
221 args.file = args.file[0]
|
|
|
222 if not args.file.lower().endswith('.pdf'):
|
|
10
|
223 sys.stderr.write("{0}: input file must end with .pdf (case insensitive)\n".format(MYNAME))
|
|
|
224 sys.exit(2)
|
|
0
|
225
|
|
|
226 # Default the set of layers to delete, if necessary
|
|
|
227 if args.layers is None:
|
|
|
228 args.layers = get_default_layers(args.file)
|
|
|
229
|
|
10
|
230 # Default the output file, if necessary
|
|
|
231 if args.output is None:
|
|
|
232 args.output = os.path.splitext(args.file)[0] + ".tiff"
|
|
|
233 elif os.path.splitext(args.output)[1].lower() not in set([".tif", ".tiff"]):
|
|
|
234 sys.stderr.write("{0}: output file must end with .tif or .tiff (case insensitive)\n".format(MYNAME))
|
|
|
235 sys.exit(2)
|
|
|
236
|
|
0
|
237 with TemporaryDirectory() as td:
|
|
|
238 # Get scratch file name. This goes under TMPDIR; if the default temporary
|
|
|
239 # area is too small, set that environment variable accordingly!
|
|
|
240 tf = os.path.join(td, "temp.tiff")
|
|
|
241
|
|
|
242 # Render. The result will have undesired margins.
|
|
|
243 sys.stdout.write("Rendering (may take a while)...\n")
|
|
|
244 sys.stdout.flush()
|
|
|
245 proc = subprocess.Popen(
|
|
|
246 gdalcmd("gdal_translate", "-q", args.file, tf, "-of", "GTiff", "--config", "GDAL_PDF_LAYERS", args.layers, "--config", "GDAL_PDF_DPI", str(args.resolution)),
|
|
|
247 stdout=subprocess.DEVNULL, stderr=subprocess.PIPE, env=ENVIRON, encoding=ENCODING)
|
|
|
248 drainout(proc.stderr)
|
|
|
249 waitfor(proc)
|
|
|
250 sys.stdout.write("Done!\n")
|
|
|
251
|
|
|
252 # Determine crop marquee.
|
|
|
253 with Image.open(tf) as im:
|
|
|
254 minx = find_edge(im, Orientation.VERTICAL, Direction.ASCENDING)
|
|
|
255 maxx = find_edge(im, Orientation.VERTICAL, Direction.DESCENDING)
|
|
|
256 miny = find_edge(im, Orientation.HORIZONTAL, Direction.ASCENDING)
|
|
|
257 maxy = find_edge(im, Orientation.HORIZONTAL, Direction.DESCENDING)
|
|
|
258
|
|
|
259 # Crop.
|
|
|
260 proc = subprocess.Popen(
|
|
10
|
261 gdalcmd("gdal_translate", "-q", tf, args.output, "-of", "GTiff", "-srcwin", str(minx), str(miny), str(maxx-minx+1), str(maxy-miny+1)),
|
|
0
|
262 stdout=subprocess.DEVNULL, stderr=sys.stderr)
|
|
|
263 waitfor(proc)
|
|
|
264
|
|
|
265 # AMF...
|
|
|
266 sys.exit(estat)
|
