← 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

  • 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
1
 
#!/bin/bash
2
 
 
3
 
############################################################################
4
 
#
5
 
# MODULE:       r.tileset for GRASS 6
6
 
# AUTHOR(S):    Cedric Shock
7
 
# PURPOSE:      To produce tilings of regions in other projections.
8
 
# COPYRIGHT:    (C) 2006 by Cedric Shoc
9
 
#
10
 
#               This program is free software under the GNU General Public
11
 
#               License (>=v2). Read the file COPYING that comes with GRASS
12
 
#               for details.
13
 
#
14
 
#############################################################################
15
 
 
16
 
# Requires:
17
 
#  bc
18
 
#  cs2cs
19
 
#  grep
20
 
#  sed
21
 
 
22
 
# Bugs:
23
 
#  Does not know about meridians in projections. However, unlike the usual
24
 
#   hack used to find meridians, this code is perfectly happy with arbitrary
25
 
#   rotations and flips
26
 
#  The following are planned to be fixed in a future version, if it turns out
27
 
#   to be necessary for someone:
28
 
#   Does not generate optimal tilings. This means that between an appropriate
29
 
#    projection for a region and latitude longitude projection, in the
30
 
#    degenerate case, it may create tiles demanding up to twice the necessary
31
 
#    information. Requesting data from cylindrical projections near their poles
32
 
#    results in divergence. You really don't want to use source data from
33
 
#    someone who put it in a cylindrical projection near a pole, do you?
34
 
#   Not generating "optimal" tilings has another side effect; the sanity
35
 
#    of the destination region will not carry over to generating tiles of
36
 
#    realistic aspect ratio. This might be a problem for some WMS servers
37
 
#    presenting data in a highly inappropriate projection. Do you really
38
 
#    want their data?
39
 
 
40
 
# Usage example
41
 
#  r.tileset -w sourceproj="+init=epsg:4326" maxcols=4096 maxrows=4096
42
 
 
43
 
# Changes:
44
 
# Version 2:
45
 
#   Added destination projection options
46
 
#   Added extra cells for overlap
47
 
 
48
 
#%Module
49
 
#% description: Produces tilings of the source projection for use in the destination region and projection.
50
 
#% keywords: raster, tiling
51
 
#%End
52
 
#%flag
53
 
#% key: g
54
 
#% description: Produces shell script output
55
 
#%end
56
 
#%flag
57
 
#% key: w
58
 
#% description: Produces web map server query string output
59
 
#%end
60
 
#%option
61
 
#% key: region
62
 
#% type: string
63
 
#% description: Name of region to use instead of current region for bounds and resolution
64
 
#%end
65
 
#%option
66
 
#% key: sourceproj
67
 
#% type: string
68
 
#% description: Source projection
69
 
#% required : yes
70
 
#%end
71
 
#%option
72
 
#% key: sourcescale
73
 
#% type: string
74
 
#% description: Conversion factor from units to meters in source projection
75
 
#% answer : 1
76
 
#%end
77
 
#%option
78
 
#% key: destproj
79
 
#% type: string
80
 
#% description: Destination projection, defaults to this location's projection
81
 
#% required : no
82
 
#%end
83
 
#%option
84
 
#% key: destscale
85
 
#% type: string
86
 
#% description: Conversion factor from units to meters in source projection
87
 
#% required : no
88
 
#%end
89
 
#%option
90
 
#% key: maxcols
91
 
#% type: integer
92
 
#% description: Maximum number of columns for a tile in the source projection
93
 
#% answer: 1024
94
 
#%end
95
 
#%option
96
 
#% key: maxrows
97
 
#% type: integer
98
 
#% description: Maximum number of rows for a tile in the source projection
99
 
#% answer: 1024
100
 
#%end
101
 
#%option
102
 
#% key: overlap
103
 
#% type: integer
104
 
#% description: Number of cells tiles should overlap in each direction
105
 
#% answer: 0
106
 
#%end
107
 
#%option
108
 
#% key: fs
109
 
#% type: string
110
 
#% description: Output field separator
111
 
#% answer:|
112
 
#%end
113
 
#%option
114
 
#% key: v
115
 
#% type: integer
116
 
#% description: Verbosity level
117
 
#% answer: 0
118
 
#%end
119
 
 
120
 
if  [ -z "$GISBASE" ] ; then
121
 
    echo "You must be in GRASS GIS to run this program." 1>&2
122
 
    exit 1
123
 
fi
124
 
 
125
 
if [ "$1" != "@ARGS_PARSED@" ] ; then
126
 
  exec g.parser "$0" "$@"
127
 
fi
128
 
 
129
 
g.message -d "[r.tileset]"
130
 
 
131
 
# Program locations:
132
 
BC="bc"
133
 
BCARGS="-l"
134
 
CS2CS="cs2cs"
135
 
 
136
 
 
137
 
# check if we have bc
138
 
if [ ! -x "`which $BC`" ] ; then
139
 
    g.message -e "$BC is required, please install it first"
140
 
    exit 1
141
 
fi
142
 
 
143
 
# check if we have cs2cs
144
 
if [ ! -x "`which $CS2CS`" ] ; then
145
 
    g.message -e "$CS2CS is required, please install it first"
146
 
    exit 1
147
 
fi
148
 
 
149
 
 
150
 
# Default IFS
151
 
defaultIFS=$IFS
152
 
 
153
 
# Data structures used in this program:
154
 
# A bounding box:
155
 
# 0 -> left, 1-> bottom, 2->right, 3-> top
156
 
# A border:
157
 
# An array of points indexed by 0 for "x" and 4 for "y" + by number 0, 1, 2, and 3
158
 
# A reprojector [0] is name of source projection, [1] is name of destination
159
 
# A projection - [0] is proj.4 text, [1] is scale
160
 
 
161
 
#####################
162
 
# name:     message
163
 
# purpose:  displays messages to the user
164
 
# usage: message level text
165
 
 
166
 
message () {
167
 
        if [ $1 -lt $GIS_OPT_V ] ; then
168
 
                shift
169
 
                echo "$@"  
170
 
        fi
171
 
}
172
 
 
173
 
#########################
174
 
# name: lookup
175
 
# purpose: lookup values in passed by reference arrays
176
 
 
177
 
Lookup_Mac() {
178
 
# Assignment Command Statement Builder
179
 
# "$1" is Source name, "$3" is Destination name
180
 
    echo "eval $3"'=${'"$1[$2]}"
181
 
}
182
 
 
183
 
# Lookup_Mac
184
 
 
185
 
declare -f LookupP               # Function "Pointer"
186
 
LookupP=Lookup_Mac               # Statement Builder
187
 
 
188
 
lookup () {
189
 
        $($LookupP $1 $2 $3)
190
 
}
191
 
 
192
 
#####################3
193
 
# name: calculate
194
 
# purpose: perform calculations
195
 
# usage: varname "expr"
196
 
 
197
 
calculate() {
198
 
        message 3 "$2"
199
 
        c_tmp=`echo "$2" | $BC $BCARGS`
200
 
        eval $1=$c_tmp
201
 
}
202
 
 
203
 
#########################
204
 
# name: makeProjector
205
 
# purpose: bundle up reprojection information
206
 
#          and provide enough abstraction to use pipes
207
 
# usage: makeReprojector source source-scale destination destination-scale name
208
 
 
209
 
makeProjector() {
210
 
        eval $5[0]=$1
211
 
        eval $5[1]=$2
212
 
        eval $5[2]=$3
213
 
        eval $5[3]=$4
214
 
}
215
 
 
216
 
#######################################################################
217
 
# name:     project
218
 
# purpose:  projects x1 y1 using projector
219
 
# usage: project -12 36 x2 y2 source-to-dest
220
 
project() {
221
 
        # echo "$5"
222
 
        # Pick apart the projector
223
 
        lookup $5 0 p_source_proj_p
224
 
        lookup $5 1 p_source_units
225
 
        lookup $5 2 p_dest_proj_p
226
 
        lookup $5 3 p_dest_units
227
 
 
228
 
        eval "p_source_proj=\$$p_source_proj_p"
229
 
        eval "p_dest_proj=\$$p_dest_proj_p"
230
 
 
231
 
        calculate p_x1 "$1 * $p_source_units"
232
 
        calculate p_y1 "$2 * $p_source_units"
233
 
 
234
 
        message 3 "echo \"$p_x1 $p_y1\" | $CS2CS -f \"%.8f\" $p_source_proj +to $p_dest_proj"
235
 
 
236
 
        answer=`echo "$p_x1 $p_y1" | eval $CS2CS -f "%.8f" $(echo $p_source_proj) +to $(echo $p_dest_proj)`
237
 
 
238
 
        if [ $? -ne 0 ] ; then
239
 
            g.message -e message="Problem running $CS2CS. <$answer>"
240
 
            exit 1
241
 
        fi
242
 
 
243
 
        message 3 "$answer"
244
 
        # do not quote $answer in the following line
245
 
        eval `echo $answer | sed -e "s/^/p_x2=/"  -e "s/ /;p_y2=/"  -e "s/ /;p_z2=/"`
246
 
 
247
 
        calculate p_x2 "$p_x2 / $p_dest_units"
248
 
        calculate p_y2 "$p_y2 / $p_dest_units"
249
 
        
250
 
        eval $3=$p_x2
251
 
        eval $4=$p_y2
252
 
}
253
 
 
254
 
min() {
255
 
        lookup "#$1" @ min_n
256
 
        lookup $1 0 min_a
257
 
        for ((min_i=0;min_i<min_n;min_i+=1)) ; do
258
 
                lookup $1 $min_i min_b
259
 
                calculate min_a "($min_b < $min_a) * $min_b + ($min_a <= $min_b) * $min_a" 
260
 
        done
261
 
        eval "$2=$min_a"
262
 
}
263
 
 
264
 
max() {
265
 
        lookup "#$1" @ max_n
266
 
        lookup $1 0 max_a
267
 
        for ((max_i=0;max_i<max_n;max_i+=1)) ; do
268
 
                lookup $1 $max_i max_b
269
 
                calculate max_a "($max_b > $max_a) * $max_b + ($max_a >= $max_b) * $max_a" 
270
 
        done
271
 
        eval "$2=$max_a"
272
 
}
273
 
 
274
 
########################
275
 
# name:    bboxToPoints
276
 
# purpose: make points that are the corners of a bounding box
277
 
 
278
 
 
279
 
bboxToPoints() {
280
 
        lookup $1 0 $2[0]
281
 
        lookup $1 1 $2[1]
282
 
        lookup $1 0 $2[2]
283
 
        lookup $1 3 $2[3]
284
 
        lookup $1 2 $2[4]
285
 
        lookup $1 3 $2[5]
286
 
        lookup $1 2 $2[6]
287
 
        lookup $1 1 $2[7]
288
 
}
289
 
 
290
 
pointsToBbox() {
291
 
        lookup $1 0 ptb_xs[0]
292
 
        lookup $1 2 ptb_xs[1]
293
 
        lookup $1 4 ptb_xs[2]
294
 
        lookup $1 6 ptb_xs[3]
295
 
        lookup $1 1 ptb_ys[0]
296
 
        lookup $1 3 ptb_ys[1]
297
 
        lookup $1 5 ptb_ys[2]
298
 
        lookup $1 7 ptb_ys[3]
299
 
        min ptb_xs $2[0] 
300
 
        min ptb_ys $2[1]
301
 
        max ptb_xs $2[2]
302
 
        max ptb_ys $2[3]
303
 
}
304
 
 
305
 
#######################################################################
306
 
# name:     projectPoints
307
 
# purpose:  projects a list of points
308
 
#
309
 
projectPoints() {
310
 
        lookup "#$1" @ pp_max
311
 
        for ((pp_i=0;pp_i<pp_max;pp_i+=2)) ; do
312
 
                calculate pp_i2 "$pp_i + 1"
313
 
                lookup $1 $pp_i pp_x1
314
 
                lookup $1 $pp_i2 pp_y1
315
 
                # echo "project $pp_x1 $pp_y1 $2[$pp_i] $2[$pp_i2] $3" 1>&2
316
 
                project $pp_x1 $pp_y1 $2[$pp_i] $2[$pp_i2] $3
317
 
        done
318
 
}
319
 
 
320
 
######################
321
 
# name:     Side Length
322
 
# purpose:  Find the length of sides of a set of points from one to the next
323
 
# usage: sideLengths points xmetric ymetric answer
324
 
 
325
 
sideLengths() {
326
 
        lookup "#$1" @ sl_max
327
 
        for ((sl_i=0;sl_i<sl_max;sl_i+=2)) ; do
328
 
                calculate sl_i2 "$sl_i + 1"
329
 
                calculate sl_j "scale=0; ($sl_i + 2) % $sl_max"
330
 
                calculate sl_j2 "$sl_j + 1"
331
 
                calculate sl_k "scale=0; $sl_i / 2"
332
 
                # echo "$sl_i, $sl_i2, $sl_j, $sl_j2"
333
 
                lookup $1 $sl_i sl_x1
334
 
                lookup $1 $sl_i2 sl_y1
335
 
                lookup $1 $sl_j sl_x2
336
 
                lookup $1 $sl_j2 sl_y2
337
 
                sl_x="(($sl_x2 - $sl_x1) * $2)"
338
 
                sl_y="(($sl_y2 - $sl_y1) * $3)"
339
 
                # echo "$sl_x, $sl_y"
340
 
                calculate sl_d "sqrt ( $sl_x * $sl_x + $sl_y * $sl_y )"
341
 
                # echo "$sl_i ($sl_k): $sl_d"
342
 
                eval "$4[$sl_k]=$sl_d"
343
 
        done
344
 
}
345
 
 
346
 
 
347
 
#######################################################################
348
 
# name:     bboxesIntersect
349
 
# purpose:  determine if two bounding boxes intersect
350
 
# usage: bboxesIntersect box1 box2 returnvar
351
 
 
352
 
bboxesIntersect() {
353
 
        lookup $1 0 bi_A1[0]
354
 
        lookup $1 1 bi_A1[1]
355
 
        lookup $1 2 bi_A2[0]
356
 
        lookup $1 3 bi_A2[1]
357
 
        lookup $2 0 bi_B1[0]
358
 
        lookup $2 1 bi_B1[1]
359
 
        lookup $2 2 bi_B2[0]
360
 
        lookup $2 3 bi_B2[1]
361
 
 
362
 
        # Fails at meridians
363
 
        # I have no way to have general knowledge of meridians
364
 
        # This means rewrite in C I imagine
365
 
        for bi_i in 0 1 ; do
366
 
                calculate IN[$bi_i] \
367
 
                "(${bi_A1[$bi_i]} <= ${bi_B1[$bi_i]} && ${bi_A2[$bi_i]} >= ${bi_B2[$bi_i]}) || \
368
 
                 (${bi_A1[$bi_i]} >= ${bi_B1[$bi_i]} && ${bi_A2[$bi_i]} <= ${bi_B2[$bi_i]}) || \
369
 
                 (${bi_A1[$bi_i]} >= ${bi_B1[$bi_i]} && ${bi_A1[$bi_i]} <= ${bi_B2[$bi_i]}) || \
370
 
                 (${bi_A2[$bi_i]} >= ${bi_B1[$bi_i]} && ${bi_A2[$bi_i]} <= ${bi_B2[$bi_i]})"
371
 
        done
372
 
 
373
 
        calculate $3 "${IN[0]} && ${IN[1]}"
374
 
}
375
 
 
376
 
### Fetch destination projection and region
377
 
SOURCE_PROJ="$GIS_OPT_SOURCEPROJ"
378
 
SOURCE_SCALE="$GIS_OPT_SOURCESCALE"
379
 
 
380
 
# Take into account those extra pixels we'll be a addin'
381
 
#MAX_COLS=$GIS_OPT_MAXCOLS
382
 
#MAX_ROWS=$GIS_OPT_MAXROWS
383
 
OVERLAP="$GIS_OPT_OVERLAP"
384
 
calculate MAX_COLS "$GIS_OPT_MAXCOLS - $OVERLAP"
385
 
calculate MAX_ROWS "$GIS_OPT_MAXROWS - $OVERLAP"
386
 
 
387
 
g.message "Getting destination projection"
388
 
 
389
 
if [ -z "$GIS_OPT_DESTPROJ" ] ; then
390
 
 
391
 
    DEST_PROJ=`g.proj -j | (
392
 
        DEST_PROJ=
393
 
        while read line ; do
394
 
            DEST_PROJ="$DEST_PROJ '$line'"
395
 
        done
396
 
        echo "$DEST_PROJ"
397
 
    )`
398
 
else
399
 
    DEST_PROJ=`echo $GIS_OPT_DESTPROJ`
400
 
fi
401
 
 
402
 
        g.message "Getting projection scale"
403
 
 
404
 
if [ -z "$GIS_OPT_DESTSCALE" ] ; then
405
 
        eval `g.proj -p | grep '^meters' | sed -e "s/:/=/" -e "s/ //g"`
406
 
        DEST_SCALE=$meters;
407
 
else
408
 
        DEST_SCALE=$GIS_OPT_DESTSCALE
409
 
fi
410
 
 
411
 
# Strip newlines from both the source and destination projections
412
 
DEST_PROJ=`echo $DEST_PROJ`
413
 
SOURCE_PROJ=`echo $SOURCE_PROJ`
414
 
 
415
 
# Set up the projections
416
 
makeProjector SOURCE_PROJ $SOURCE_SCALE DEST_PROJ $DEST_SCALE SOURCE_TO_DEST
417
 
makeProjector DEST_PROJ $DEST_SCALE SOURCE_PROJ $SOURCE_SCALE DEST_TO_SOURCE
418
 
 
419
 
g.message "Getting destination region"
420
 
if [ -z $GIS_OPT_REGION ] ; then
421
 
        eval `g.region -g`
422
 
else
423
 
        eval `g.region region=$GIS_OPT_REGION -g`
424
 
fi
425
 
 
426
 
 
427
 
DEST_BBOX[0]=$w
428
 
DEST_BBOX[1]=$s
429
 
DEST_BBOX[2]=$e
430
 
DEST_BBOX[3]=$n
431
 
DEST_NSRES=$nsres
432
 
DEST_EWRES=$ewres
433
 
calculate DEST_ROWS "($n - $s) / $nsres"
434
 
calculate DEST_COLS "($e - $w) / $ewres"
435
 
 
436
 
message 1 "Rows: $DEST_ROWS, Cols: $DEST_COLS"
437
 
 
438
 
### Project the destination region into the source:
439
 
g.message "Projecting destination region into source"
440
 
 
441
 
bboxToPoints DEST_BBOX DEST_BBOX_POINTS
442
 
 
443
 
message 1 "${DEST_BBOX_POINTS[@]}"
444
 
 
445
 
# echo "${DEST_TO_SOURCE[@]}"
446
 
# echo "${DEST_TO_SOURCE[0]}"
447
 
# echo "${DEST_TO_SOURCE[1]}"
448
 
# echo "${DEST_TO_SOURCE[2]}"
449
 
# echo "${DEST_TO_SOURCE[3]}"
450
 
 
451
 
# exit 1
452
 
 
453
 
projectPoints DEST_BBOX_POINTS DEST_BBOX_SOURCE_POINTS DEST_TO_SOURCE 
454
 
 
455
 
message 1 "${DEST_BBOX_SOURCE_POINTS[@]}"
456
 
 
457
 
pointsToBbox DEST_BBOX_SOURCE_POINTS SOURCE_BBOX
458
 
 
459
 
message 1 "${SOURCE_BBOX[@]}"
460
 
 
461
 
g.message "Projecting source bounding box into destination"
462
 
 
463
 
bboxToPoints SOURCE_BBOX SOURCE_BBOX_POINTS
464
 
 
465
 
message 1 "${SOURCE_BBOX_POINTS[@]}"
466
 
 
467
 
projectPoints SOURCE_BBOX_POINTS SOURCE_BBOX_DEST_POINTS SOURCE_TO_DEST 
468
 
 
469
 
message 1 "${SOURCE_BBOX_DEST_POINTS[@]}"
470
 
 
471
 
calculate X_METRIC "1 / $DEST_EWRES"
472
 
calculate Y_METRIC "1 / $DEST_NSRES"
473
 
 
474
 
g.message "Computing length of sides of source bounding box"
475
 
 
476
 
sideLengths  SOURCE_BBOX_DEST_POINTS $X_METRIC $Y_METRIC SOURCE_BBOX_DEST_LENGTHS
477
 
 
478
 
message 1 "${SOURCE_BBOX_DEST_LENGTHS[@]}"
479
 
 
480
 
# Find the skewedness of the two directions.
481
 
# Define it to be greater than one
482
 
# In the direction (x or y) in which the world is least skewed (ie north south in lat long)
483
 
# Divide the world into strips. These strips are as big as possible contrained by max_
484
 
# In the other direction do the same thing.
485
 
# Theres some recomputation of the size of the world that's got to come in here somewhere.
486
 
 
487
 
# For now, however, we are going to go ahead and request more data than is necessary.
488
 
# For small regions far from the critical areas of projections this makes very little difference
489
 
# in the amount of data gotten.
490
 
# We can make this efficient for big regions or regions near critical points later.
491
 
 
492
 
HORIZONTALS=( ${SOURCE_BBOX_DEST_LENGTHS[1]} ${SOURCE_BBOX_DEST_LENGTHS[3]} )
493
 
VERTICALS=( ${SOURCE_BBOX_DEST_LENGTHS[0]} ${SOURCE_BBOX_DEST_LENGTHS[2]} )
494
 
 
495
 
max HORIZONTALS BIGGER[0]
496
 
max VERTICALS BIGGER[1]
497
 
 
498
 
MAX[0]=$MAX_COLS
499
 
MAX[1]=$MAX_ROWS
500
 
 
501
 
# Compute the number and size of tiles to use in each direction
502
 
# I'm making fairly even sized tiles
503
 
# They differer from each other in height and width only by one cell
504
 
# I'm going to make the numbers all simpler and add this extra cell to
505
 
# every tile.
506
 
 
507
 
g.message "Computing tiling"
508
 
 
509
 
for i in 0 1 ; do
510
 
        # Make these into integers.
511
 
        # Round up
512
 
        calculate BIGGER[$i] "scale=0; (${BIGGER[$i]} + 1) / 1"
513
 
        calculate TILES[$i] "scale=0; (${BIGGER[$i]} / ${MAX[$i]}) + 1"
514
 
        # BC truncates and doesn't round, which is perfect here:
515
 
        calculate TILE_BASE_SIZE[$i] "scale=0; (${BIGGER[$i]} / ${TILES[$i]})"
516
 
        calculate TILES_EXTRA_1[$i] "scale=0; (${BIGGER[$i]} % ${TILES[$i]})"
517
 
        # This is adding the extra pixel (remainder) to all of the tiles:
518
 
        calculate TILE_SIZE[$i] "scale=0; ${TILE_BASE_SIZE[$i]} + (${TILES_EXTRA_1[$i]} > 0)" 
519
 
        calculate TILESET_SIZE[$i] "scale=0; ${TILE_SIZE[$i]} * ${TILES[$i]}"
520
 
        # Add overlap to tiles (doesn't effect tileset_size
521
 
        calculate TILE_SIZE_OVERLAP[$i] "${TILE_SIZE[$i]} + $OVERLAP"
522
 
done;
523
 
 
524
 
g.message "There will be ${TILES[0]} by ${TILES[1]} tiles each ${TILE_SIZE[0]} by ${TILE_SIZE[1]} cells."
525
 
 
526
 
ximax=${TILES[0]}
527
 
yimax=${TILES[1]}
528
 
 
529
 
MIN_X=${SOURCE_BBOX[0]}
530
 
MIN_Y=${SOURCE_BBOX[1]}
531
 
MAX_X=${SOURCE_BBOX[2]}
532
 
MAX_Y=${SOURCE_BBOX[3]}
533
 
SPAN_X="($MAX_X - $MIN_X)"
534
 
SPAN_Y="($MAX_Y - $MIN_Y)"
535
 
 
536
 
 
537
 
for ((xi=0;xi<ximax;xi+=1)); do
538
 
        calculate TILE_BBOX[0] "$MIN_X + ( $xi * ${TILE_SIZE[0]} / ${TILESET_SIZE[0]} ) * $SPAN_X"
539
 
        calculate TILE_BBOX[2] "$MIN_X + ( ($xi + 1) * ${TILE_SIZE_OVERLAP[0]} / ${TILESET_SIZE[0]} ) * $SPAN_X"
540
 
        for ((yi=0;yi<yimax;yi+=1)); do
541
 
                calculate TILE_BBOX[1] "$MIN_Y + ( $yi * ${TILE_SIZE[1]} / ${TILESET_SIZE[1]} ) * $SPAN_Y"
542
 
                calculate TILE_BBOX[3] "$MIN_Y + ( ($yi + 1) * ${TILE_SIZE_OVERLAP[1]} / ${TILESET_SIZE[1]} ) * $SPAN_Y"
543
 
                bboxToPoints TILE_BBOX TILE_BBOX_POINTS
544
 
                projectPoints TILE_BBOX_POINTS TILE_DEST_BBOX_POINTS SOURCE_TO_DEST
545
 
                pointsToBbox TILE_DEST_BBOX_POINTS TILE_DEST_BBOX
546
 
                bboxesIntersect TILE_DEST_BBOX DEST_BBOX intersect
547
 
                if [ $intersect ] ; then
548
 
                        if [ $GIS_FLAG_W -eq 1 ] ; then
549
 
                                echo "bbox=${TILE_BBOX[0]},${TILE_BBOX[1]},${TILE_BBOX[2]},${TILE_BBOX[3]}&width=${TILE_SIZE_OVERLAP[0]}&height=${TILE_SIZE_OVERLAP[1]}"
550
 
                        elif [ $GIS_FLAG_G -eq 1 ] ; then
551
 
                                echo "w=${TILE_BBOX[0]};s=${TILE_BBOX[1]};e=${TILE_BBOX[2]};n=${TILE_BBOX[3]};cols=${TILE_SIZE_OVERLAP[0]};rows=${TILE_SIZE_OVERLAP[1]};"
552
 
                        else
553
 
                                echo "${TILE_BBOX[0]}$GIS_OPT_FS${TILE_BBOX[1]}$GIS_OPT_FS${TILE_BBOX[2]}$GIS_OPT_FS${TILE_BBOX[3]}$GIS_OPT_FS${TILE_SIZE_OVERLAP[0]}$GIS_OPT_FS${TILE_SIZE_OVERLAP[1]}"
554
 
                        fi
555
 
                fi 
556
 
        done
557
 
done
558