source: subversion/applications/utils/export/osm2pgsql-intarray/expire-tiles.c @ 28719

Last change on this file since 28719 was 19148, checked in by jonb, 10 years ago

osm2pgsql: Apply multipolygon patch from Twain with a few changes.

File size: 13.5 KB
Line 
1/*
2 * Dirty tile list generation
3 *
4 * Steve Hill <steve@nexusuk.org>
5 *
6 * Please refer to the OpenPisteMap expire_tiles.py script for a demonstration
7 * of how to make use of the output:
8 * https://subversion.nexusuk.org/trac/browser/openpistemap/trunk/scripts/expire_tiles.py
9 */
10
11#include <libpq-fe.h>
12#include <math.h>
13#include <stdlib.h>
14#include <string.h>
15#include <errno.h>
16#include "expire-tiles.h"
17#include "output.h"
18#include "pgsql.h"
19#include "build_geometry.h"
20#include "reprojection.h"
21
22#define EARTH_CIRCUMFERENCE             40075016.68
23#define HALF_EARTH_CIRCUMFERENCE        (EARTH_CIRCUMFERENCE / 2)
24#define TILE_EXPIRY_LEEWAY              0.5             // How many tiles worth of space to leave either side of a changed feature
25#define EXPIRE_TILES_MAX_BBOX           30000           // Maximum width or height of a bounding box (metres)
26
27struct tile {
28        int             complete[2][2]; // Flags
29        struct tile *   subtiles[2][2];
30};
31
32int map_width; /* not "static" since used in reprojection.c! */
33static double                           tile_width;
34static const struct output_options *    Options;
35static struct tile *                    dirty = NULL;
36static int                              outcount;
37
38/*
39 * We store the dirty tiles in an in-memory tree during runtime
40 * and dump them out to a file at the end.  This allows us to easilly drop
41 * duplicate tiles from the output.
42 *
43 * This data structure consists of a node, representing a tile at zoom level 0,
44 * which contains 4 pointers to nodes representing each of the child tiles at
45 * zoom level 1, and so on down the the zoom level specified in
46 * Options->expire_tiles_zoom.
47 *
48 * The memory allowed to this structure is not capped, but daily deltas
49 * generally produce a few hundred thousand expired tiles at zoom level 17,
50 * which are easilly accommodated.
51 */
52
53static int calc_complete(struct tile * tile) {
54        int     c;
55
56        c = tile->complete[0][0];
57        c += tile->complete[0][1];
58        c += tile->complete[1][0];
59        c += tile->complete[1][1];
60        return c;
61}
62
63static void destroy_tree(struct tile * tree) {
64        if (! tree) return;
65        if (tree->subtiles[0][0]) destroy_tree(tree->subtiles[0][0]);
66        if (tree->subtiles[0][1]) destroy_tree(tree->subtiles[0][1]);
67        if (tree->subtiles[1][0]) destroy_tree(tree->subtiles[1][0]);
68        if (tree->subtiles[1][1]) destroy_tree(tree->subtiles[1][1]);
69        free(tree);
70}
71
72/*
73 * Mark a tile as dirty.
74 * Returns the number of subtiles which have all their children marked as dirty.
75 */
76static int _mark_tile(struct tile ** tree, int x, int y, int zoom, int this_zoom) {
77        int     zoom_diff = zoom - this_zoom;
78        int     rel_x;
79        int     rel_y;
80        int     complete;
81
82        if (! *tree) *tree = calloc(1, sizeof(**tree));
83        zoom_diff = (zoom - this_zoom) - 1;
84        rel_x = (x >> zoom_diff) & 1;
85        rel_y = (y >> zoom_diff) & 1;
86        if (! (*tree)->complete[rel_x][rel_y]) {
87                if (zoom_diff <= 0) {
88                        (*tree)->complete[rel_x][rel_y] = 1;
89                } else {
90                        complete = _mark_tile(&((*tree)->subtiles[rel_x][rel_y]), x, y, zoom, this_zoom + 1);
91                        if (complete >= 4) {
92                                (*tree)->complete[rel_x][rel_y] = 1;
93                                // We can destroy the subtree to save memory now all the children are dirty
94                                destroy_tree((*tree)->subtiles[rel_x][rel_y]);
95                                (*tree)->subtiles[rel_x][rel_y] = NULL;
96                        }
97                }
98        }
99        return calc_complete(*tree);
100}
101
102/*
103 * Mark a tile as dirty.
104 * Returns the number of subtiles which have all their children marked as dirty.
105 */
106static int mark_tile(struct tile ** tree_head, int x, int y, int zoom) {
107        return _mark_tile(tree_head, x, y, zoom, 0);
108}
109
110static void output_dirty_tile(FILE * outfile, int x, int y, int zoom, int min_zoom) {
111        int     y_min;
112        int     x_iter;
113        int     y_iter;
114        int     x_max;
115        int     y_max;
116        int     out_zoom;
117        int     zoom_diff;
118
119        if (zoom > min_zoom) out_zoom = zoom;
120        else out_zoom = min_zoom;
121        zoom_diff = out_zoom - zoom;
122        y_min = y << zoom_diff;
123        x_max = (x + 1) << zoom_diff;
124        y_max = (y + 1) << zoom_diff;
125        for (x_iter = x << zoom_diff; x_iter < x_max; x_iter++) {
126                for (y_iter = y_min; y_iter < y_max; y_iter++) {
127                        outcount++;
128                        if ((outcount <= 1) || (! (outcount % 1000))) {
129                                fprintf(stderr, "\rWriting dirty tile list (%iK)", outcount / 1000);
130                                fflush(stderr);
131                        }
132                        fprintf(outfile, "%i/%i/%i\n", out_zoom, x_iter, y_iter);
133                }
134        }
135}
136
137static void _output_and_destroy_tree(FILE * outfile, struct tile * tree, int x, int y, int this_zoom, int min_zoom) {
138        int     sub_x = x << 1;
139        int     sub_y = y << 1;
140        FILE *  ofile;
141
142        if (! tree) return;
143
144        ofile = outfile;
145        if ((tree->complete[0][0]) && outfile) {
146                output_dirty_tile(outfile, sub_x + 0, sub_y + 0, this_zoom + 1, min_zoom);
147                ofile = NULL;
148        }
149        if (tree->subtiles[0][0]) _output_and_destroy_tree(ofile, tree->subtiles[0][0], sub_x + 0, sub_y + 0, this_zoom + 1, min_zoom);
150
151        ofile = outfile;
152        if ((tree->complete[0][1]) && outfile) {
153                output_dirty_tile(outfile, sub_x + 0, sub_y + 1, this_zoom + 1, min_zoom);
154                ofile = NULL;
155        }
156        if (tree->subtiles[0][1]) _output_and_destroy_tree(ofile, tree->subtiles[0][1], sub_x + 0, sub_y + 1, this_zoom + 1, min_zoom);
157
158        ofile = outfile;
159        if ((tree->complete[1][0]) && outfile) {
160                output_dirty_tile(outfile, sub_x + 1, sub_y + 0, this_zoom + 1, min_zoom);
161                ofile = NULL;
162        }
163        if (tree->subtiles[1][0]) _output_and_destroy_tree(ofile, tree->subtiles[1][0], sub_x + 1, sub_y + 0, this_zoom + 1, min_zoom);
164
165        ofile = outfile;
166        if ((tree->complete[1][1]) && outfile) {
167                output_dirty_tile(outfile, sub_x + 1, sub_y + 1, this_zoom + 1, min_zoom);
168                ofile = NULL;
169        }
170        if (tree->subtiles[1][1]) _output_and_destroy_tree(ofile, tree->subtiles[1][1], sub_x + 1, sub_y + 1, this_zoom + 1, min_zoom);
171
172        free(tree);
173}
174
175static void output_and_destroy_tree(FILE * outfile, struct tile * tree) {
176        _output_and_destroy_tree(outfile, tree, 0, 0, 0, Options->expire_tiles_zoom_min);
177}
178
179void expire_tiles_stop(void) {
180        FILE *  outfile;
181
182        if (Options->expire_tiles_zoom < 0) return;
183        outcount = 0;
184        if ((outfile = fopen(Options->expire_tiles_filename, "a"))) {
185            output_and_destroy_tree(outfile, dirty);
186            fclose(outfile);
187        } else {
188        fprintf(stderr, "Failed to open expired tiles file (%s).  Tile expiry list will not be written!\n", strerror(errno));
189    }
190        dirty = NULL;
191}
192
193void expire_tiles_init(const struct output_options *options) {
194        Options = options;
195        if (Options->expire_tiles_zoom < 0) return;
196        map_width = 1 << Options->expire_tiles_zoom;
197        tile_width = EARTH_CIRCUMFERENCE / map_width;
198}
199
200static void expire_tile(int x, int y) {
201        mark_tile(&dirty, x, y, Options->expire_tiles_zoom);
202}
203
204static int normalise_tile_x_coord(int x) {
205        x %= map_width;
206        if (x < 0) x = (map_width - x) + 1;
207        return x;
208}
209
210/*
211 * Expire tiles that a line crosses
212 */
213static void expire_tiles_from_line(double lon_a, double lat_a, double lon_b, double lat_b) {
214        double  tile_x_a;
215        double  tile_y_a;
216        double  tile_x_b;
217        double  tile_y_b;
218        double  temp;
219        double  x1;
220        double  y1;
221        double  x2;
222        double  y2;
223        double  hyp_len;
224        double  x_len;
225        double  y_len;
226        double  x_step;
227        double  y_step;
228        double  step;
229        double  next_step;
230        int     x;
231        int     y;
232        int     norm_x;
233
234    coords_to_tile(&tile_x_a, &tile_y_a, lon_a, lat_a);
235    coords_to_tile(&tile_x_b, &tile_y_b, lon_b, lat_b);
236
237        if (tile_x_a > tile_x_b) {
238                // We always want the line to go from left to right - swap the ends if it doesn't
239                temp = tile_x_b;
240                tile_x_b = tile_x_a;
241                tile_x_a = temp;
242                temp = tile_y_b;
243                tile_y_b = tile_y_a;
244                tile_y_a = temp;
245        }
246
247        x_len = tile_x_b - tile_x_a;
248        if (x_len > map_width / 2) {
249                // If the line is wider than half the map, assume it
250                // crosses the international date line.
251                // These coordinates get normalised again later
252                tile_x_a += map_width;
253                temp = tile_x_b;
254                tile_x_b = tile_x_a;
255                tile_x_a = temp;
256                temp = tile_y_b;
257                tile_y_b = tile_y_a;
258                tile_y_a = temp;
259        }
260        y_len = tile_y_b - tile_y_a;
261        hyp_len = sqrt(pow(x_len, 2) + pow(y_len, 2));  // Pythagoras
262        x_step = x_len / hyp_len;
263        y_step = y_len / hyp_len;
264//      fprintf(stderr, "Expire from line (%f,%f),(%f,%f) [%f,%f],[%f,%f] %fx%f hyp_len = %f\n", lon_a, lat_a, lon_b, lat_b, tile_x_a, tile_y_a, tile_x_b, tile_y_b, x_len, y_len, hyp_len);
265       
266        for (step = 0; step <= hyp_len; step ++) {
267                // Interpolate points 1 tile width apart
268                next_step = step + 1;
269                if (next_step > hyp_len) next_step = hyp_len;
270                x1 = tile_x_a + ((double)step * x_step);
271                y1 = tile_y_a + ((double)step * y_step);
272                x2 = tile_x_a + ((double)next_step * x_step);
273                y2 = tile_y_a + ((double)next_step * y_step);
274               
275//              printf("Expire from subline (%f,%f),(%f,%f)\n", x1, y1, x2, y2);
276                // The line (x1,y1),(x2,y2) is up to 1 tile width long
277                // x1 will always be <= x2
278                // We could be smart and figure out the exact tiles intersected,
279                // but for simplicity, treat the coordinates as a bounding box
280                // and expire everything within that box.
281                if (y1 > y2) {
282                        temp = y2;
283                        y2 = y1;
284                        y1 = temp;
285                }
286                for (x = x1 - TILE_EXPIRY_LEEWAY; x <= x2 + TILE_EXPIRY_LEEWAY; x ++) {
287                        norm_x =  normalise_tile_x_coord(x);
288                        for (y = y1 - TILE_EXPIRY_LEEWAY; y <= y2 + TILE_EXPIRY_LEEWAY; y ++) {
289                                expire_tile(norm_x, y);
290                        }
291                }
292        }
293}
294
295/*
296 * Expire tiles within a bounding box
297 */
298int expire_tiles_from_bbox(double min_lon, double min_lat, double max_lon, double max_lat) {
299        double          width;
300        double          height;
301        int             min_tile_x;
302        int             min_tile_y;
303        int             max_tile_x;
304        int             max_tile_y;
305        int             iterator_x;
306        int             iterator_y;
307        int             norm_x;
308        int             ret;
309    double  tmp_x;
310    double  tmp_y;
311
312        if (Options->expire_tiles_zoom < 0) return 0;
313
314        width = max_lon - min_lon;
315        height = max_lat - min_lat;
316        if (width > HALF_EARTH_CIRCUMFERENCE + 1) {
317                // Over half the planet's width within the bounding box - assume the
318                // box crosses the international date line and split it into two boxes
319                ret = expire_tiles_from_bbox(-HALF_EARTH_CIRCUMFERENCE, min_lat, min_lon, max_lat);
320                ret += expire_tiles_from_bbox(max_lon, min_lat, HALF_EARTH_CIRCUMFERENCE, max_lat);
321                return ret;
322        }
323
324        if (width > EXPIRE_TILES_MAX_BBOX) return -1;
325        if (height > EXPIRE_TILES_MAX_BBOX) return -1;
326
327//      printf("Expire from bbox (%f,%f)-(%f,%f) %fx%f\n", min_lon, min_lat, min_lon, min_lat, width, height);
328
329        // Convert the box's Mercator coordinates into tile coordinates
330    coords_to_tile(&tmp_x, &tmp_y, min_lon, max_lat);
331    min_tile_x = tmp_x - TILE_EXPIRY_LEEWAY;
332    min_tile_y = tmp_y - TILE_EXPIRY_LEEWAY;
333    coords_to_tile(&tmp_x, &tmp_y, max_lon, min_lat);
334    max_tile_x = tmp_x + TILE_EXPIRY_LEEWAY;
335    max_tile_y = tmp_y + TILE_EXPIRY_LEEWAY;
336        if (min_tile_x < 0) min_tile_x = 0;
337        if (min_tile_y < 0) min_tile_y = 0;
338        if (max_tile_x > map_width) max_tile_x = map_width;
339        if (max_tile_y > map_width) max_tile_y = map_width;
340//      printf("BBOX: (%f %f) - (%f %f) [%i %i] - [%i %i]\n", min_lon, min_lat, max_lon, max_lat, min_tile_x, min_tile_y, max_tile_x, max_tile_y);
341        for (iterator_x = min_tile_x; iterator_x <= max_tile_x; iterator_x ++) {
342                norm_x =  normalise_tile_x_coord(iterator_x);
343                for (iterator_y = min_tile_y; iterator_y <= max_tile_y; iterator_y ++) {
344                        expire_tile(norm_x, iterator_y);
345                }
346        }
347        return 0;
348}
349
350void expire_tiles_from_nodes_line(struct osmNode * nodes, int count) {
351        int     i;
352        double  last_lat;
353        double  last_lon;
354
355        if (Options->expire_tiles_zoom < 0) return;
356//      fprintf(stderr, "Expire from nodes_line (%i)\n", count);
357        if (count < 1) return;
358        last_lat = nodes[0].lat;
359        last_lon = nodes[0].lon;
360        if (count < 2) {
361                expire_tiles_from_bbox(last_lon, last_lat, last_lon, last_lat);
362                return;
363        }
364        for (i = 1; i < count; i ++) {
365                expire_tiles_from_line(last_lon, last_lat, nodes[i].lon, nodes[i].lat);
366                last_lat = nodes[i].lat;
367                last_lon = nodes[i].lon;
368        }
369}
370
371/*
372 * Calculate a bounding box from a list of nodes and expire all tiles within it
373 */
374void expire_tiles_from_nodes_poly(struct osmNode * nodes, int count, int osm_id) {
375        int     i;
376        int     got_coords = 0;
377        double  min_lon = 0.0;
378        double  min_lat = 0.0;
379        double  max_lon = 0.0;
380        double  max_lat = 0.0;
381       
382        if (Options->expire_tiles_zoom < 0) return;
383//      printf("Expire from nodes_poly (%i)\n", count);
384        for (i = 0; i < count; i++) {
385                if ((! got_coords) || (nodes[i].lon < min_lon)) min_lon = nodes[i].lon;
386                if ((! got_coords) || (nodes[i].lat < min_lat)) min_lat = nodes[i].lat;
387                if ((! got_coords) || (nodes[i].lon > max_lon)) max_lon = nodes[i].lon;
388                if ((! got_coords) || (nodes[i].lat > max_lat)) max_lat = nodes[i].lat;
389                got_coords = 1;
390        }
391        if (got_coords) {
392                if (expire_tiles_from_bbox(min_lon, min_lat, max_lon, max_lat)) {
393                        // Bounding box too big - just expire tiles on the line
394                        fprintf(stderr, "\rLarge polygon (%.0f x %.0f metres, OSM ID %i) - only expiring perimeter\n", max_lon - min_lon, max_lat - min_lat, osm_id);
395                        expire_tiles_from_nodes_line(nodes, count);
396                }
397        }
398}
399
400static void expire_tiles_from_xnodes_poly(struct osmNode ** xnodes, int * xcount, int osm_id) {
401        int     i;
402
403        for (i = 0; xnodes[i]; i++) expire_tiles_from_nodes_poly(xnodes[i], xcount[i], osm_id);
404}
405
406static void expire_tiles_from_xnodes_line(struct osmNode ** xnodes, int * xcount) {
407        int     i;
408
409        for (i = 0; xnodes[i]; i++) expire_tiles_from_nodes_line(xnodes[i], xcount[i]);
410}
411
412void expire_tiles_from_wkt(const char * wkt, int osm_id) {
413        struct osmNode **       xnodes;
414        int *                   xcount;
415        int                     polygon;
416        int                     i;
417
418        if (Options->expire_tiles_zoom < 0) return;
419        if (! parse_wkt(wkt, &xnodes, &xcount, &polygon)) {
420                if (polygon) expire_tiles_from_xnodes_poly(xnodes, xcount, osm_id);
421                else expire_tiles_from_xnodes_line(xnodes, xcount);
422                for (i = 0; xnodes[i]; i++) free(xnodes[i]);
423                free(xnodes);
424                free(xcount);
425        }
426}
427
428void expire_tiles_from_db(PGconn * sql_conn, int osm_id) {
429        PGresult *      res;
430        char            tmp[16];
431        char const *    paramValues[1];
432        int             tuple;
433        char *          wkt;
434
435        if (Options->expire_tiles_zoom < 0) return;
436        snprintf(tmp, sizeof(tmp), "%d", osm_id);
437        paramValues[0] = tmp;
438        res = pgsql_execPrepared(sql_conn, "get_way", 1, paramValues, PGRES_TUPLES_OK);
439        for (tuple = 0; tuple < PQntuples(res); tuple++) {
440                wkt = PQgetvalue(res, tuple, 0);
441                expire_tiles_from_wkt(wkt, osm_id);
442        }
443        PQclear(res);
444}
445
446
Note: See TracBrowser for help on using the repository browser.