Tools/Prep/ShapeFile/shape.cxx

   1 // shape.cxx -- shape/gpc utils
   2 //
   3 // Written by Curtis Olson, started February 1999.
   4 //
   5 // Copyright (C) 1999  Curtis L. Olson  - curt@flightgear.org
   6 //
   7 // This program is free software; you can redistribute it and/or modify
   8 // it under the terms of the GNU General Public License as published by
   9 // the Free Software Foundation; either version 2 of the License, or
  10 // (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful,
  13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 // GNU General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20 //
  21 // $Id$
  22
  23
  24 #include <Include/compiler.h>
  25
  26 #include STL_STRING
  27
  28 #include <Bucket/newbucket.hxx>
  29 #include <Debug/logstream.hxx>
  30
  31 #include <Polygon/index.hxx>
  32 #include <Polygon/names.hxx>
  33 #include "shape.hxx"
  34
  35
  36 #define FG_MAX_VERTICES 100000
  37 static gpc_vertex_list v_list;
  38
  39
  40 class point2d {
  41 public:
  42     double x, y;
  43 };
  44
  45
  46 static void clip_and_write_poly( string root, long int p_index, AreaType area,
  47                                  FGBucket b, gpc_polygon *shape ) {
  48     point2d c, min, max;
  49     c.x = b.get_center_lon();
  50     c.y = b.get_center_lat();
  51     double span = bucket_span(c.y);
  52     gpc_polygon base, result;
  53     char tile_name[256], poly_index[256];
  54
  55     // calculate bucket dimensions
  56     if ( (c.y >= -89.0) && (c.y < 89.0) ) {
  57         min.x = c.x - span / 2.0;
  58         max.x = c.x + span / 2.0;
  59         min.y = c.y - FG_HALF_BUCKET_SPAN;
  60         max.y = c.y + FG_HALF_BUCKET_SPAN;
  61     } else if ( c.y < -89.0) {
  62         min.x = -90.0;
  63         max.x = -89.0;
  64         min.y = -180.0;
  65         max.y = 180.0;
  66     } else if ( c.y >= 89.0) {
  67         min.x = 89.0;
  68         max.x = 90.0;
  69         min.y = -180.0;
  70         max.y = 180.0;
  71     } else {
  72         FG_LOG ( FG_GENERAL, FG_ALERT,
  73                  "Out of range latitude in clip_and_write_poly() = " << c.y );
  74     }
  75
  76     FG_LOG( FG_GENERAL, FG_INFO, "  (" << min.x << "," << min.y << ") ("
  77             << max.x << "," << max.y << ")" );
  78
  79     // set up clipping tile
  80     v_list.vertex[0].x = min.x;
  81     v_list.vertex[0].y = min.y;
  82
  83     v_list.vertex[1].x = max.x;
  84     v_list.vertex[1].y = min.y;
  85
  86     v_list.vertex[2].x = max.x;
  87     v_list.vertex[2].y = max.y;
  88
  89     v_list.vertex[3].x = min.x;
  90     v_list.vertex[3].y = max.y;
  91
  92     v_list.num_vertices = 4;
  93
  94     base.num_contours = 0;
  95     base.contour = NULL;
  96     base.hole = NULL;
  97     gpc_add_contour( &base, &v_list, 0 );
  98
  99     // FG_LOG( FG_GENERAL, FG_DEBUG, "base = 4 vertices" );
 100
 101     /*
 102     FILE *bfp= fopen("base", "w");
 103     gpc_write_polygon(bfp, &base);
 104     fclose(bfp);
 105     */
 106
 107     gpc_polygon_clip(GPC_INT, &base, shape, &result);
 108
 109     if ( result.num_contours > 0 ) {
 110         long int t_index = b.gen_index();
 111         string path = root + "/Scenery/" + b.gen_base_path();
 112         string command = "mkdir -p " + path;
 113         system( command.c_str() );
 114
 115         sprintf( tile_name, "%ld", t_index );
 116         string polyfile = path + "/" + tile_name;
 117
 118         sprintf( poly_index, "%ld", p_index );
 119         polyfile += ".";
 120         polyfile += poly_index;
 121
 122         string poly_type = get_area_name( area );
 123         if ( poly_type == "Unknown" ) {
 124             cout << "unknown area type in clip_and_write_poly()!" << endl;
 125             exit(-1);
 126         }
 127
 128         FILE *rfp= fopen( polyfile.c_str(), "w" );
 129         fprintf( rfp, "%s\n", poly_type.c_str() );
 130         gpc_write_polygon( rfp, 1, &result );
 131         fclose( rfp );
 132
 133         // only free result if it is not empty
 134         gpc_free_polygon(&result);
 135     }
 136
 137     gpc_free_polygon(&base);
 138 }
 139
 140
 141 // Initialize structure we use to create polygons for the gpc library
 142 bool shape_utils_init() {
 143     v_list.num_vertices = 0;
 144     v_list.vertex = new gpc_vertex[FG_MAX_VERTICES];;
 145
 146     return true;
 147 }
 148
 149
 150 // initialize a gpc_polygon
 151 void init_shape(gpc_polygon *shape) {
 152     shape->num_contours = 0;
 153     shape->contour = NULL;
 154     shape->hole = NULL;
 155 }
 156
 157
 158 // make a gpc_polygon, first contour is outline, remaining contours
 159 // are holes
 160 void add_to_shape(int count, double *coords, gpc_polygon *shape) {
 161
 162     for ( int i = 0; i < count; i++ ) {
 163         v_list.vertex[i].x = coords[i*2+0];
 164         v_list.vertex[i].y = coords[i*2+1];
 165     }
 166
 167     v_list.num_vertices = count;
 168
 169     if ( shape->num_contours == 0 ) {
 170         // outline
 171         gpc_add_contour( shape, &v_list, 0 );
 172     } else {
 173         // hole
 174         gpc_add_contour( shape, &v_list, 1 );
 175     }
 176 }
 177
 178
 179 // process shape (write polygon to all intersecting tiles)
 180 void process_shape(string path, AreaType area, gpc_polygon *master_gpc_shape) {
 181     point2d min, max;
 182     long int index;
 183     int i, j;
 184
 185     min.x = min.y = 200.0;
 186     max.x = max.y = -200.0;
 187
 188     // find min/max of polygon
 189     for ( i = 0; i < master_gpc_shape->num_contours; i++ ) {
 190         for ( j = 0; j < master_gpc_shape->contour[i].num_vertices; j++ ) {
 191             double x = master_gpc_shape->contour[i].vertex[j].x;
 192             double y = master_gpc_shape->contour[i].vertex[j].y;
 193
 194             if ( x < min.x ) { min.x = x; }
 195             if ( y < min.y ) { min.y = y; }
 196             if ( x > max.x ) { max.x = x; }
 197             if ( y > max.y ) { max.y = y; }
 198         }
 199     }
 200
 201     /*
 202     FILE *sfp= fopen("shape", "w");
 203     gpc_write_polygon(sfp, master_gpc_shape);
 204     fclose(sfp);
 205     exit(-1);
 206     */
 207
 208     // get next polygon index
 209     index = poly_index_next();
 210
 211     FG_LOG( FG_GENERAL, FG_INFO, "  min = " << min.x << "," << min.y
 212             << " max = " << max.x << "," << max.y );
 213
 214     // find buckets for min, and max points of convex hull.
 215     // note to self: self, you should think about checking for
 216     // polygons that span the date line
 217     FGBucket b_min(min.x, min.y);
 218     FGBucket b_max(max.x, max.y);
 219     FG_LOG( FG_GENERAL, FG_INFO, "  Bucket min = " << b_min );
 220     FG_LOG( FG_GENERAL, FG_INFO, "  Bucket max = " << b_max );
 221
 222     if ( b_min == b_max ) {
 223         clip_and_write_poly( path, index, area, b_min, master_gpc_shape );
 224     } else {
 225         FGBucket b_cur;
 226         int dx, dy, i, j;
 227
 228         fgBucketDiff(b_min, b_max, &dx, &dy);
 229         FG_LOG( FG_GENERAL, FG_INFO,
 230                 "  polygon spans tile boundaries" );
 231         FG_LOG( FG_GENERAL, FG_INFO, "  dx = " << dx
 232                 << "  dy = " << dy );
 233
 234         if ( (dx > 200) || (dy > 200) ) {
 235             FG_LOG( FG_GENERAL, FG_ALERT,
 236                     "somethings really wrong!!!!" );
 237             exit(-1);
 238         }
 239
 240         for ( j = 0; j <= dy; j++ ) {
 241             // for performance reasons, we'll clip out just this
 242             // horizontal row, and clip all the tiles in this row
 243             // against the smaller shape
 244
 245             FG_LOG ( FG_GENERAL, FG_INFO,
 246                      "Generating clip row " << j << " of " << dy );
 247
 248             FGBucket b_clip = fgBucketOffset(min.x, min.y, 0, j);
 249             gpc_polygon row, clip_row;
 250             point2d c, clip_max, clip_min;
 251             c.x = b_clip.get_center_lon();
 252             c.y = b_clip.get_center_lat();
 253
 254             // calculate bucket clip_min.y and clip_max.y
 255             if ( (c.y >= -89.0) && (c.y < 89.0) ) {
 256                 clip_min.y = c.y - FG_HALF_BUCKET_SPAN;
 257                 clip_max.y = c.y + FG_HALF_BUCKET_SPAN;
 258             } else if ( c.y < -89.0) {
 259                 clip_min.y = -90.0;
 260                 clip_max.y = -89.0;
 261             } else if ( c.y >= 89.0) {
 262                 clip_min.y = 89.0;
 263                 clip_max.y = 90.0;
 264             } else {
 265                 FG_LOG ( FG_GENERAL, FG_ALERT,
 266                          "Out of range latitude in clip_and_write_poly() = "
 267                          << c.y );
 268             }
 269             clip_min.x = -180.0;
 270             clip_max.x = 180.0;
 271
 272             // set up clipping tile
 273             v_list.vertex[0].x = clip_min.x;
 274             v_list.vertex[0].y = clip_min.y;
 275
 276             v_list.vertex[1].x = clip_max.x;
 277             v_list.vertex[1].y = clip_min.y;
 278
 279             v_list.vertex[2].x = clip_max.x;
 280             v_list.vertex[2].y = clip_max.y;
 281
 282             v_list.vertex[3].x = clip_min.x;
 283             v_list.vertex[3].y = clip_max.y;
 284
 285             v_list.num_vertices = 4;
 286
 287             row.num_contours = 0;
 288             row.contour = NULL;
 289             row.hole = NULL;
 290             gpc_add_contour( &row, &v_list, 0 );
 291
 292             clip_row.num_contours = 0;
 293             clip_row.contour = NULL;
 294             clip_row.hole = NULL;
 295
 296             gpc_polygon_clip(GPC_INT, &row, master_gpc_shape, &clip_row);
 297
 298             /* FILE *sfp = fopen("master_shape", "w");
 299             gpc_write_polygon(sfp, 0, master_gpc_shape);
 300             fclose(sfp);
 301             sfp = fopen("clip_row", "w");
 302             gpc_write_polygon(sfp, 0, &clip_row);
 303             fclose(sfp); */
 304
 305             for ( i = 0; i <= dx; i++ ) {
 306                 b_cur = fgBucketOffset(min.x, min.y, i, j);
 307                 clip_and_write_poly( path, index, area, b_cur, &clip_row );
 308             }
 309             gpc_free_polygon(&row);
 310             gpc_free_polygon(&clip_row);
 311         }
 312         // string answer; cin >> answer;
 313     }
 314 }
 315
 316
 317 // free a gpc_polygon
 318 void free_shape(gpc_polygon *shape) {
 319     gpc_free_polygon(shape);
 320 }
 321
 322