Mercurial > cgi-bin > hgweb.cgi > TopoTiler

comparison make-tile @ 0:c20b5314774f

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