← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to scripts/r.tileset/r.tileset.py

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/r.tileset/r.tileset.py
 
1
#!/usr/bin/env python
 
2
 
 
3
############################################################################
 
4
#
 
5
# MODULE:       r.tileset
 
6
#
 
7
# AUTHOR(S):    Cedric Shock
 
8
#               Updated for GRASS7 by Martin Landa, 2009
 
9
#
 
10
# PURPOSE:      To produce tilings of regions in other projections.
 
11
#
 
12
# COPYRIGHT:    (C) 2006-2009 by Cedric Shoc, Martin Landa, and GRASS development team
 
13
#
 
14
#               This program is free software under the GNU General
 
15
#               Public License (>=v2). Read the file COPYING that
 
16
#               comes with GRASS for details.
 
17
#
 
18
#############################################################################
 
19
 
 
20
# Bugs:
 
21
#  Does not know about meridians in projections. However, unlike the usual
 
22
#   hack used to find meridians, this code is perfectly happy with arbitrary
 
23
#   rotations and flips
 
24
#  The following are planned to be fixed in a future version, if it turns out
 
25
#   to be necessary for someone:
 
26
#   Does not generate optimal tilings. This means that between an appropriate
 
27
#    projection for a region and latitude longitude projection, in the
 
28
#    degenerate case, it may create tiles demanding up to twice the necessary
 
29
#    information. Requesting data from cylindrical projections near their poles
 
30
#    results in divergence. You really don't want to use source data from
 
31
#    someone who put it in a cylindrical projection near a pole, do you?
 
32
#   Not generating "optimal" tilings has another side effect; the sanity
 
33
#    of the destination region will not carry over to generating tiles of
 
34
#    realistic aspect ratio. This might be a problem for some WMS servers
 
35
#    presenting data in a highly inappropriate projection. Do you really
 
36
#    want their data?
 
37
 
 
38
#%module
 
39
#% description: Produces tilings of the source projection for use in the destination region and projection.
 
40
#% keyword: raster
 
41
#% keyword: tiling
 
42
#%end
 
43
#%flag
 
44
#% key: g
 
45
#% description: Produces shell script output
 
46
#%end
 
47
#%flag
 
48
#% key: w
 
49
#% description: Produces web map server query string output
 
50
#%end
 
51
#%option
 
52
#% key: region
 
53
#% type: string
 
54
#% description: Name of region to use instead of current region for bounds and resolution
 
55
#%end
 
56
#%option
 
57
#% key: sourceproj
 
58
#% type: string
 
59
#% description: Source projection
 
60
#% required : yes
 
61
#%end
 
62
#%option
 
63
#% key: sourcescale
 
64
#% type: string
 
65
#% description: Conversion factor from units to meters in source projection
 
66
#% answer : 1
 
67
#%end
 
68
#%option
 
69
#% key: destproj
 
70
#% type: string
 
71
#% description: Destination projection, defaults to this location's projection
 
72
#% required : no
 
73
#%end
 
74
#%option
 
75
#% key: destscale
 
76
#% type: string
 
77
#% description: Conversion factor from units to meters in source projection
 
78
#% required : no
 
79
#%end
 
80
#%option
 
81
#% key: maxcols
 
82
#% type: integer
 
83
#% description: Maximum number of columns for a tile in the source projection
 
84
#% answer: 1024
 
85
#%end
 
86
#%option
 
87
#% key: maxrows
 
88
#% type: integer
 
89
#% description: Maximum number of rows for a tile in the source projection
 
90
#% answer: 1024
 
91
#%end
 
92
#%option
 
93
#% key: overlap
 
94
#% type: integer
 
95
#% description: Number of cells tiles should overlap in each direction
 
96
#% answer: 0
 
97
#%end
 
98
#%option G_OPT_F_SEP
 
99
#% description: Output field separator
 
100
#%end
 
101
 
 
102
# Data structures used in this program:
 
103
# A bounding box:
 
104
# 0 -> left, 1-> bottom, 2->right, 3-> top
 
105
# A border:
 
106
# An array of points indexed by 0 for "x" and 4 for "y" + by number 0, 1, 2, and 3
 
107
# A reprojector [0] is name of source projection, [1] is name of destination
 
108
# A projection - [0] is proj.4 text, [1] is scale
 
109
 
 
110
import sys
 
111
import tempfile
 
112
import math
 
113
 
 
114
from grass.script.utils import separator
 
115
from grass.script import core as grass
 
116
from grass.exceptions import CalledModuleError
 
117
 
 
118
 
 
119
def bboxToPoints(bbox):
 
120
    """Make points that are the corners of a bounding box"""
 
121
    points = []
 
122
    points.append((bbox['w'], bbox['s']))
 
123
    points.append((bbox['w'], bbox['n']))
 
124
    points.append((bbox['e'], bbox['n']))
 
125
    points.append((bbox['e'], bbox['s']))
 
126
 
 
127
    return points
 
128
 
 
129
 
 
130
def pointsToBbox(points):
 
131
    bbox = {}
 
132
    min_x = min_y = max_x = max_y = None
 
133
    for point in points:
 
134
        if not min_x:
 
135
            min_x = max_x = point[0]
 
136
        if not min_y:
 
137
            min_y = max_y = point[1]
 
138
 
 
139
        if min_x > point[0]:
 
140
            min_x = point[0]
 
141
        if max_x < point[0]:
 
142
            max_x = point[0]
 
143
        if min_y > point[1]:
 
144
            min_y = point[1]
 
145
        if max_y < point[1]:
 
146
            max_y = point[1]
 
147
 
 
148
    bbox['n'] = max_y
 
149
    bbox['s'] = min_y
 
150
    bbox['w'] = min_x
 
151
    bbox['e'] = max_x
 
152
 
 
153
    return bbox
 
154
 
 
155
 
 
156
def project(file, source, dest):
 
157
    """Projects point (x, y) using projector"""
 
158
    errors = 0
 
159
    points = []
 
160
    try:
 
161
        ret = grass.read_command('m.proj',
 
162
                                 quiet=True,
 
163
                                 flags='d',
 
164
                                 proj_in=source['proj'],
 
165
                                 proj_out=dest['proj'],
 
166
                                 sep=';',
 
167
                                 input=file)
 
168
    except CalledModuleError:
 
169
        grass.fatal(cs2cs + ' failed')
 
170
 
 
171
    if not ret:
 
172
        grass.fatal(cs2cs + ' failed')
 
173
 
 
174
    for line in ret.splitlines():
 
175
        if "*" in line:
 
176
            errors += 1
 
177
        else:
 
178
            p_x2, p_y2, p_z2 = map(float, line.split(';'))
 
179
            points.append((p_x2 / dest['scale'], p_y2 / dest['scale']))
 
180
 
 
181
    return points, errors
 
182
 
 
183
 
 
184
def projectPoints(points, source, dest):
 
185
    """Projects a list of points"""
 
186
    dest_points = []
 
187
 
 
188
    input = tempfile.NamedTemporaryFile(mode="wt")
 
189
    for point in points:
 
190
        input.file.write('%f;%f\n' % (point[0] * source['scale'],
 
191
                                      point[1] * source['scale']))
 
192
    input.file.flush()
 
193
 
 
194
    dest_points, errors = project(input.name, source, dest)
 
195
 
 
196
    return dest_points, errors
 
197
 
 
198
 
 
199
def sideLengths(points, xmetric, ymetric):
 
200
    """Find the length of sides of a set of points from one to the next"""
 
201
    ret = []
 
202
    for i in range(len(points)):
 
203
        x1, y1 = points[i]
 
204
        j = i + 1
 
205
        if j >= len(points):
 
206
            j = 0
 
207
        sl_x = (points[j][0] - points[i][0]) * xmetric
 
208
        sl_y = (points[j][1] - points[i][1]) * ymetric
 
209
        sl_d = math.sqrt(sl_x * sl_x + sl_y * sl_y)
 
210
        ret.append(sl_d)
 
211
 
 
212
    return {'x': (ret[1], ret[3]), 'y': (ret[0], ret[2])}
 
213
 
 
214
 
 
215
def bboxesIntersect(bbox_1, bbox_2):
 
216
    """Determine if two bounding boxes intersect"""
 
217
    bi_a1 = (bbox_1['w'], bbox_1['s'])
 
218
    bi_a2 = (bbox_1['e'], bbox_1['n'])
 
219
    bi_b1 = (bbox_2['w'], bbox_2['s'])
 
220
    bi_b2 = (bbox_2['e'], bbox_2['n'])
 
221
    cin = [False, False]
 
222
    for i in (0, 1):
 
223
        if (bi_a1[i] <= bi_b1[i] and bi_a2[i] >= bi_b1[i]) or \
 
224
           (bi_a1[i] <= bi_b1[i] and bi_a2[i] >= bi_b2[i]) or \
 
225
           (bi_b1[i] <= bi_a1[i] and bi_b2[i] >= bi_a1[i]) or \
 
226
           (bi_b1[i] <= bi_a1[i] and bi_b2[i] >= bi_a2[i]):
 
227
            cin[i] = True
 
228
 
 
229
    if cin[0] and cin[1]:
 
230
        return True
 
231
 
 
232
    return False
 
233
 
 
234
 
 
235
def main():
 
236
    # Take into account those extra pixels we'll be a addin'
 
237
    max_cols = int(options['maxcols']) - int(options['overlap'])
 
238
    max_rows = int(options['maxrows']) - int(options['overlap'])
 
239
 
 
240
    if max_cols == 0:
 
241
        grass.fatal(_("It is not possibile to set 'maxcols=%s' and "
 
242
                      "'overlap=%s'. Please set maxcols>overlap" %
 
243
                      (options['maxcols'], options['overlap'])))
 
244
    elif max_rows == 0:
 
245
        grass.fatal(_("It is not possibile to set 'maxrows=%s' and "
 
246
                      "'overlap=%s'. Please set maxrows>overlap" %
 
247
                      (options['maxrows'], options['overlap'])))
 
248
    # destination projection
 
249
    if not options['destproj']:
 
250
        dest_proj = grass.read_command('g.proj',
 
251
                                       quiet=True,
 
252
                                       flags='jf').rstrip('\n')
 
253
        if not dest_proj:
 
254
            grass.fatal(_('g.proj failed'))
 
255
    else:
 
256
        dest_proj = options['destproj']
 
257
    grass.debug("Getting destination projection -> '%s'" % dest_proj)
 
258
 
 
259
    # projection scale
 
260
    if not options['destscale']:
 
261
        ret = grass.parse_command('g.proj',
 
262
                                  quiet=True,
 
263
                                  flags='j')
 
264
        if not ret:
 
265
            grass.fatal(_('g.proj failed'))
 
266
 
 
267
        if '+to_meter' in ret:
 
268
            dest_scale = ret['+to_meter'].strip()
 
269
        else:
 
270
            grass.warning(_("Scale (%s) not found, assuming '1'") % '+to_meter')
 
271
            dest_scale = '1'
 
272
    else:
 
273
        dest_scale = options['destscale']
 
274
    grass.debug('Getting destination projection scale -> %s' % dest_scale)
 
275
 
 
276
    # set up the projections
 
277
    srs_source = {'proj': options['sourceproj'],
 
278
                  'scale': float(options['sourcescale'])}
 
279
    srs_dest = {'proj': dest_proj, 'scale': float(dest_scale)}
 
280
 
 
281
    if options['region']:
 
282
        grass.run_command('g.region',
 
283
                          quiet=True,
 
284
                          region=options['region'])
 
285
    dest_bbox = grass.region()
 
286
    grass.debug('Getting destination region')
 
287
 
 
288
    # output field separator
 
289
    fs = separator(options['separator'])
 
290
 
 
291
    # project the destination region into the source:
 
292
    grass.verbose('Projecting destination region into source...')
 
293
    dest_bbox_points = bboxToPoints(dest_bbox)
 
294
 
 
295
    dest_bbox_source_points, errors_dest = projectPoints(dest_bbox_points,
 
296
                                                         source=srs_dest,
 
297
                                                         dest=srs_source)
 
298
 
 
299
    if len(dest_bbox_source_points) == 0:
 
300
        grass.fatal(_("There are no tiles available. Probably the output "
 
301
                      "projection system it is not compatible with the "
 
302
                      "projection of the current location"))
 
303
 
 
304
    source_bbox = pointsToBbox(dest_bbox_source_points)
 
305
 
 
306
    grass.verbose('Projecting source bounding box into destination...')
 
307
 
 
308
    source_bbox_points = bboxToPoints(source_bbox)
 
309
 
 
310
    source_bbox_dest_points, errors_source = projectPoints(source_bbox_points,
 
311
                                                            source=srs_source,
 
312
                                                            dest=srs_dest)
 
313
 
 
314
    x_metric = 1 / dest_bbox['ewres']
 
315
    y_metric = 1 / dest_bbox['nsres']
 
316
 
 
317
    grass.verbose('Computing length of sides of source bounding box...')
 
318
 
 
319
    source_bbox_dest_lengths = sideLengths(source_bbox_dest_points,
 
320
                                           x_metric, y_metric)
 
321
 
 
322
    # Find the skewedness of the two directions.
 
323
    # Define it to be greater than one
 
324
    # In the direction (x or y) in which the world is least skewed (ie north south in lat long)
 
325
    # Divide the world into strips. These strips are as big as possible contrained by max_
 
326
    # In the other direction do the same thing.
 
327
    # Theres some recomputation of the size of the world that's got to come in here somewhere.
 
328
 
 
329
    # For now, however, we are going to go ahead and request more data than is necessary.
 
330
    # For small regions far from the critical areas of projections this makes very little difference
 
331
    # in the amount of data gotten.
 
332
    # We can make this efficient for big regions or regions near critical points later.
 
333
 
 
334
    bigger = []
 
335
    bigger.append(max(source_bbox_dest_lengths['x']))
 
336
    bigger.append(max(source_bbox_dest_lengths['y']))
 
337
    maxdim = (max_cols, max_rows)
 
338
 
 
339
    # Compute the number and size of tiles to use in each direction
 
340
    # I'm making fairly even sized tiles
 
341
    # They differer from each other in height and width only by one cell
 
342
    # I'm going to make the numbers all simpler and add this extra cell to
 
343
    # every tile.
 
344
 
 
345
    grass.message(_('Computing tiling...'))
 
346
    tiles = [-1, -1]
 
347
    tile_base_size = [-1, -1]
 
348
    tiles_extra_1 = [-1, -1]
 
349
    tile_size = [-1, -1]
 
350
    tileset_size = [-1, -1]
 
351
    tile_size_overlap = [-1, -1]
 
352
    for i in range(len(bigger)):
 
353
        # make these into integers.
 
354
        # round up
 
355
        bigger[i] = int(bigger[i] + 1)
 
356
        tiles[i] = int((bigger[i] / maxdim[i]) + 1)
 
357
        tile_size[i] = tile_base_size[i] = int(bigger[i] / tiles[i])
 
358
        tiles_extra_1[i] = int(bigger[i] % tiles[i])
 
359
        # This is adding the extra pixel (remainder) to all of the tiles:
 
360
        if tiles_extra_1[i] > 0:
 
361
            tile_size[i] = tile_base_size[i] + 1
 
362
        tileset_size[i] = int(tile_size[i] * tiles[i])
 
363
        # Add overlap to tiles (doesn't effect tileset_size
 
364
        tile_size_overlap[i] = tile_size[i] + int(options['overlap'])
 
365
 
 
366
    grass.verbose("There will be %d by %d tiles each %d by %d cells" %
 
367
                  (tiles[0], tiles[1], tile_size[0], tile_size[1]))
 
368
 
 
369
    ximax = tiles[0]
 
370
    yimax = tiles[1]
 
371
 
 
372
    min_x = source_bbox['w']
 
373
    min_y = source_bbox['s']
 
374
    max_x = source_bbox['e']
 
375
    max_y = source_bbox['n']
 
376
    span_x = (max_x - min_x)
 
377
    span_y = (max_y - min_y)
 
378
 
 
379
    xi = 0
 
380
    tile_bbox = {'w': -1, 's': -1, 'e': -1, 'n': -1}
 
381
 
 
382
    if errors_dest > 0:
 
383
        grass.warning(_("During computation %i tiles could not be created" %
 
384
                        errors_dest))
 
385
 
 
386
    while xi < ximax:
 
387
        tile_bbox['w'] = float(min_x) + (float(xi) * float(tile_size[0]) / float(tileset_size[0])) * float(span_x)
 
388
        tile_bbox['e'] = float(min_x) + (float(xi + 1) * float(tile_size_overlap[0]) / float(tileset_size[0])) * float(span_x)
 
389
        yi = 0
 
390
        while yi < yimax:
 
391
            tile_bbox['s'] = float(min_y) + (float(yi) * float(tile_size[1]) / float(tileset_size[1])) * float(span_y)
 
392
            tile_bbox['n'] = float(min_y) + (float(yi + 1) * float(tile_size_overlap[1]) / float(tileset_size[1])) * float(span_y)
 
393
            tile_bbox_points = bboxToPoints(tile_bbox)
 
394
            tile_dest_bbox_points, errors = projectPoints(tile_bbox_points,
 
395
                                                          source=srs_source,
 
396
                                                          dest=srs_dest)
 
397
            tile_dest_bbox = pointsToBbox(tile_dest_bbox_points)
 
398
            if bboxesIntersect(tile_dest_bbox, dest_bbox):
 
399
                if flags['w']:
 
400
                    print "bbox=%s,%s,%s,%s&width=%s&height=%s" % \
 
401
                          (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'],
 
402
                           tile_bbox['n'], tile_size_overlap[0],
 
403
                           tile_size_overlap[1])
 
404
                elif flags['g']:
 
405
                    print "w=%s;s=%s;e=%s;n=%s;cols=%s;rows=%s" % \
 
406
                          (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'],
 
407
                           tile_bbox['n'], tile_size_overlap[0],
 
408
                           tile_size_overlap[1])
 
409
                else:
 
410
                    print "%s%s%s%s%s%s%s%s%s%s%s" % \
 
411
                          (tile_bbox['w'], fs, tile_bbox['s'], fs,
 
412
                           tile_bbox['e'], fs, tile_bbox['n'], fs,
 
413
                           tile_size_overlap[0], fs, tile_size_overlap[1])
 
414
            yi += 1
 
415
        xi += 1
 
416
 
 
417
if __name__ == "__main__":
 
418
    cs2cs = 'cs2cs'
 
419
    options, flags = grass.parser()
 
420
    sys.exit(main())