GIS/InstallationGuidelines/Linux: osm2pgsql_centroid.patch

File osm2pgsql_centroid.patch, 18.3 KB (added by Fran Boon, 14 years ago)

  • output.h

     
    2929  const char *expire_tiles_filename;    /* File name to output expired tiles list to */
    3030  int enable_hstore; /* add an additional hstore column with objects key/value pairs */
    3131  int enable_multi; /* Output multi-geometries intead of several simple geometries */
     32  int enable_centroid; /* Generate an additional geometry column that contains the center of the geometry */
    3233  const char** hstore_columns; /* list of columns that should be written into their own hstore column */
    3334  int n_hstore_columns; /* number of hstore columns */
    3435};

  • build_geometry.cpp

     
    2121*/
    2222
    2323#include <iostream>
     24#include <cstdio>
    2425#include <cstring>
    2526#include <cstdlib>
    2627#include <exception>
    … …  
    6061typedef std::auto_ptr<Geometry> geom_ptr;
    6162
    6263static std::vector<std::string> wkts;
     64static std::vector<std::string> wkt_centroids;
    6365static std::vector<double> areas;
    6466
    6567
    … …  
    8183            std::auto_ptr<LinearRing> shell(gf.createLinearRing(coords.release()));
    8284            geom = geom_ptr(gf.createPolygon(shell.release(), new std::vector<Geometry *>));
    8385            geom->normalize(); // Fix direction of ring
     86           
     87            // this function seems not to be used anymore so centroid functionality is not implemented here
     88            //   see get_wkt_split
    8489        } else {
    8590            if (coords->getSize() < 2)
    8691                return NULL;
    … …  
    104109    }
    105110}
    106111
    107 
    108 size_t get_wkt_split(osmNode *nodes, int count, int polygon, double split_at) {
     112size_t get_wkt_split(osmNode *nodes, int count, int polygon, double split_at, int enable_centroid) {
    109113    GeometryFactory gf;
    110114    std::auto_ptr<CoordinateSequence> coords(gf.getCoordinateSequenceFactory()->create(0, 2));
    111115    double area;
    112116    WKTWriter wktw;
    113117    size_t wkt_size = 0;
     118    char centroid[100];
    114119
    115120    try
    116121    {
    … …  
    130135            std::string wkt = wktw.write(geom.get());
    131136            wkts.push_back(wkt);
    132137            areas.push_back(area);
     138           
     139            // centroid of a closed way
     140            if(enable_centroid)
     141            {
     142                Point* p = geom->getInteriorPoint();
     143                snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)", p->getX(), p->getY());
     144                wkt_centroids.push_back(centroid);
     145            }
     146           
    133147            wkt_size++;
    134148        } else {
    135149            if (coords->getSize() < 2)
    … …  
    147161                    std::string wkt = wktw.write(geom.get());
    148162                    wkts.push_back(wkt);
    149163                    areas.push_back(0);
     164                   
     165                    // centroid of a linear way
     166                    if(enable_centroid)
     167                    {
     168                        // the interior point is always on an edge of the line
     169                        Point* p = geom->getInteriorPoint();
     170                       
     171                        // the centroid is not always part of the line
     172                        // Point* p = geom->getCentroid();
     173                       
     174                        // interpolate a point on the geometry, needs geos >= 3.2
     175                        // Point* p = geom->interpolate(0.5);
     176                       
     177                        snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)", p->getX(), p->getY());
     178                        wkt_centroids.push_back(centroid);
     179                    }
     180                   
    150181                    wkt_size++;
    151182                    distance=0;
    152183                    segment = std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0, 2));
    … …  
    182213        return result;
    183214}
    184215
     216char * get_wkt_centroid(size_t index)
     217{
     218    if( wkt_centroids.size() <= index )
     219        return NULL;
     220
     221    char *result;
     222    result = (char*) std::malloc( wkt_centroids[index].length() + 1);
     223   
     224    // At least give some idea of why we about to seg fault
     225    if (!result) std::cerr << std::endl << "Unable to allocate memory: " << (wkt_centroids[index].length() + 1) << std::endl;
     226   
     227    std::strcpy(result, wkt_centroids[index].c_str());
     228    return result;
     229}
     230
    185231double get_area(size_t index)
    186232{
    187233    return areas[index];
    … …  
    191237{
    192238   wkts.clear();
    193239   areas.clear();
     240   wkt_centroids.clear();
    194241}
    195242
    196243static int coords2nodes(CoordinateSequence * coords, struct osmNode ** nodes) {
    … …  
    291338    return 1;
    292339}
    293340
    294 size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon, int enable_multi, double split_at) {
     341size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon, int enable_multi, double split_at, int enable_centroid) {
    295342    size_t wkt_size = 0;
    296343    std::auto_ptr<std::vector<Geometry*> > lines(new std::vector<Geometry*>);
    297344    GeometryFactory gf;
    298345    geom_ptr geom;
    299 
     346    char centroid[100];
     347   
    300348    try
    301349    {
    302350        for (int c=0; xnodes[c]; c++) {
    … …  
    359407                        std::string wkt = writer.write(geom.get());
    360408                        wkts.push_back(wkt);
    361409                        areas.push_back(0);
     410                       
     411                        // centroid of a route-relation (linear)
     412                        if(enable_centroid)
     413                        {
     414                            Point* p = geom->getInteriorPoint();
     415                            snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)", p->getX(), p->getY());
     416                            wkt_centroids.push_back(centroid);
     417                        }
     418                       
    362419                        wkt_size++;
    363420                        distance=0;
    364421                        segment = std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0, 2));
    … …  
    457514                std::string text = writer.write(multipoly.get());
    458515                wkts.push_back(text);
    459516                areas.push_back(multipoly->getArea());
     517               
     518                // centroid of a multipolygon relation (area) with more
     519                // then one outer way in in multi-geometry mode
     520                if(enable_centroid)
     521                {
     522                    Point* p = multipoly->getInteriorPoint();
     523                    snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)", p->getX(), p->getY());
     524                    wkt_centroids.push_back(centroid);
     525                }
     526               
    460527                wkt_size++;
    461528            }
    462529            else
    … …  
    467534                    std::string text = writer.write(poly);
    468535                    wkts.push_back(text);
    469536                    areas.push_back(poly->getArea());
     537                   
     538                    // centroid of a multipolygon relation (area)
     539                    if(enable_centroid)
     540                    {
     541                        Point* p = poly->getInteriorPoint();
     542                        snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)", p->getX(), p->getY());
     543                        wkt_centroids.push_back(centroid);
     544                    }
     545                   
    470546                    wkt_size++;
    471547                    delete(poly);
    472548                }

  • osm2pgsql.c

     
    180180    printf("                     \tthat start with the specified string, eg --hstore-column \"name:\" will\n");
    181181    printf("                     \tproduce an extra hstore column that contains all name:xx tags\n");
    182182    printf("   -G|--multi-geometry\t\tGenerate multi-geometry features in postgresql tables.\n");
     183    printf("   -n|--centroid\t\tGenerate an additional geometry column that contains the center of the geometry.\n");
    183184    printf("   -h|--help\t\tHelp information.\n");
    184185    printf("   -v|--verbose\t\tVerbose output.\n");
    185186    printf("\n");
    … …  
    299300    int expire_tiles_zoom_min = -1;
    300301    int enable_hstore = 0;
    301302    int enable_multi = 0;
     303    int enable_centroid = 0;
    302304    const char *expire_tiles_filename = "dirty_tiles";
    303305    const char *db = "gis";
    304306    const char *username=NULL;
    … …  
    352354            {"hstore", 0, 0, 'k'},
    353355            {"hstore-column", 1, 0, 'z'},
    354356            {"multi-geometry", 0, 0, 'G'},
     357            {"centroid", 0, 0, 'n'},
    355358            {"input-reader", 1, 0, 'r'},
    356359            {"version", 0, 0, 'V'},
    357360            {0, 0, 0, 0}
    … …  
    397400                hstore_columns[n_hstore_columns-1] = optarg;
    398401                break;
    399402            case 'G': enable_multi=1; break;
     403        case 'n': enable_centroid=1; break;
    400404            case 'r': input_reader = optarg; break;
    401405            case 'h': long_usage_bool=1; break;
    402406            case 'V': exit(EXIT_SUCCESS);
    … …  
    469473    options.expire_tiles_zoom_min = expire_tiles_zoom_min;
    470474    options.expire_tiles_filename = expire_tiles_filename;
    471475    options.enable_multi = enable_multi;
     476    options.enable_centroid = enable_centroid;
    472477    options.enable_hstore = enable_hstore;
    473478    options.hstore_columns = hstore_columns;
    474479    options.n_hstore_columns = n_hstore_columns;

  • build_geometry.h

     
    3232int parse_wkt(const char * wkt, struct osmNode *** xnodes, int ** xcount, int * polygon);
    3333
    3434char *get_wkt_simple(struct osmNode *, int count, int polygon);
    35 size_t get_wkt_split(struct osmNode *, int count, int polygon, double split_at);
     35size_t get_wkt_split(struct osmNode *, int count, int polygon, double split_at, int enable_centroid);
    3636
    3737char* get_wkt(size_t index);
     38char * get_wkt_centroid(size_t index);
    3839double get_area(size_t index);
    39 size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon, int enable_multi, double split_at);
     40size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int make_polygon, int enable_multi, double split_at, int enable_centroid);
    4041void clear_wkts();
    4142
    4243#ifdef __cplusplus

  • output-pgsql.c

     
    5151    unsigned int buflen;
    5252    int copyMode;
    5353    char *columns;
     54    int hasCentroidColumn;
    5455} tables [] = {
    5556    { .name = "%s_point",   .type = "POINT"     },
    5657    { .name = "%s_line",    .type = "LINESTRING"},
    … …  
    265266    /* Return to copy mode if we dropped out */
    266267    if( !tables[table].copyMode )
    267268    {
    268         pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN", tables[table].name, tables[table].columns);
     269        if( tables[table].hasCentroidColumn == 1 )
     270            pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way,centroid) FROM STDIN", tables[table].name, tables[table].columns);
     271        else
     272            pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN", tables[table].name, tables[table].columns);
     273       
    269274        tables[table].copyMode = 1;
    270275    }
    271276    /* If the combination of old and new data is too big, flush old data */
    … …  
    631636
    632637
    633638
    634 static void write_wkts(int id, struct keyval *tags, const char *wkt, enum table_id table)
     639static void write_wkts(int id, struct keyval *tags, const char *wkt, enum table_id table, const char *wkt_centroid)
    635640{
    636641 
    637642    static char *sql;
    … …  
    674679    sprintf(sql, "SRID=%d;", SRID);
    675680    copy_to_table(table, sql);
    676681    copy_to_table(table, wkt);
     682   
     683    if(wkt_centroid && tables[table].hasCentroidColumn == 1)
     684    {
     685        sprintf(sql, "\tSRID=%d;", SRID);
     686        copy_to_table(table, sql);
     687        copy_to_table(table, wkt_centroid);
     688    }
     689   
    677690    copy_to_table(table, "\n");
    678691}
    679692
    … …  
    800813    else
    801814        split_at = 100 * 1000;
    802815
    803     wkt_size = get_wkt_split(nodes, count, polygon, split_at);
     816    wkt_size = get_wkt_split(nodes, count, polygon, split_at, Options->enable_centroid);
    804817
    805818    for (i=0;i<wkt_size;i++)
    806819    {
    807820        char *wkt = get_wkt(i);
     821        char *wkt_centroid = NULL;
    808822
     823        if(Options->enable_centroid)
     824            wkt_centroid = get_wkt_centroid(i);
     825
    809826        if (wkt && strlen(wkt)) {
    810827            /* FIXME: there should be a better way to detect polygons */
    811828            if (!strncmp(wkt, "POLYGON", strlen("POLYGON")) || !strncmp(wkt, "MULTIPOLYGON", strlen("MULTIPOLYGON"))) {
    … …  
    816833                    snprintf(tmp, sizeof(tmp), "%f", area);
    817834                    addItem(tags, "way_area", tmp, 0);
    818835                }
    819                 write_wkts(id, tags, wkt, t_poly);
     836                write_wkts(id, tags, wkt, t_poly, wkt_centroid);
    820837            } else {
    821838                expire_tiles_from_nodes_line(nodes, count);
    822                 write_wkts(id, tags, wkt, t_line);
     839                write_wkts(id, tags, wkt, t_line, wkt_centroid);
    823840                if (roads)
    824                     write_wkts(id, tags, wkt, t_roads);
     841                    write_wkts(id, tags, wkt, t_roads, wkt_centroid);
    825842            }
    826843        }
    827844        free(wkt);
     845        if(wkt_centroid) free(wkt_centroid);
    828846    }
    829847    clear_wkts();
    830848       
    … …  
    10261044    else
    10271045        split_at = 100 * 1000;
    10281046
    1029     wkt_size = build_geometry(id, xnodes, xcount, make_polygon, Options->enable_multi, split_at);
     1047    wkt_size = build_geometry(id, xnodes, xcount, make_polygon, Options->enable_multi, split_at, Options->enable_centroid);
    10301048
    10311049    if (!wkt_size) {
    10321050        resetList(&tags);
    … …  
    10371055    for (i=0;i<wkt_size;i++)
    10381056    {
    10391057        char *wkt = get_wkt(i);
     1058        char *wkt_centroid = NULL;
    10401059
     1060        if(Options->enable_centroid)
     1061            wkt_centroid = get_wkt_centroid(i);
     1062
    10411063        if (wkt && strlen(wkt)) {
    10421064            expire_tiles_from_wkt(wkt, -id);
    10431065            /* FIXME: there should be a better way to detect polygons */
    … …  
    10481070                    snprintf(tmp, sizeof(tmp), "%f", area);
    10491071                    addItem(&tags, "way_area", tmp, 0);
    10501072                }
    1051                 write_wkts(-id, &tags, wkt, t_poly);
     1073                write_wkts(-id, &tags, wkt, t_poly, wkt_centroid);
    10521074            } else {
    1053                 write_wkts(-id, &tags, wkt, t_line);
     1075                write_wkts(-id, &tags, wkt, t_line, wkt_centroid);
    10541076                if (roads)
    1055                     write_wkts(-id, &tags, wkt, t_roads);
     1077                    write_wkts(-id, &tags, wkt, t_roads, wkt_centroid);
    10561078            }
    10571079        }
    10581080        free(wkt);
    … …  
    10881110    // If we are making a boundary then also try adding any relations which form complete rings
    10891111    // The linear variants will have already been processed above
    10901112    if (make_boundary) {
    1091         wkt_size = build_geometry(id, xnodes, xcount, 1, Options->enable_multi, split_at);
     1113        wkt_size = build_geometry(id, xnodes, xcount, 1, Options->enable_multi, split_at, Options->enable_centroid);
    10921114        for (i=0;i<wkt_size;i++)
    10931115        {
    10941116            char *wkt = get_wkt(i);
     1117            char *wkt_centroid = NULL;
    10951118
     1119            if(Options->enable_centroid)
     1120                wkt_centroid = get_wkt_centroid(i);
     1121
    10961122            if (strlen(wkt)) {
    10971123                expire_tiles_from_wkt(wkt, -id);
    10981124                /* FIXME: there should be a better way to detect polygons */
    … …  
    11031129                        snprintf(tmp, sizeof(tmp), "%f", area);
    11041130                        addItem(&tags, "way_area", tmp, 0);
    11051131                    }
    1106                     write_wkts(-id, &tags, wkt, t_poly);
     1132                    write_wkts(-id, &tags, wkt, t_poly, wkt_centroid);
    11071133                }
    11081134            }
    11091135            free(wkt);
    … …  
    12101236            pgsql_exec(sql_conn, PGRES_COMMAND_OK, "%s", sql);
    12111237            pgsql_exec(sql_conn, PGRES_TUPLES_OK, "SELECT AddGeometryColumn('%s', 'way', %d, '%s', 2 );\n",
    12121238                        tables[i].name, SRID, tables[i].type );
     1239           
     1240            // if a centroid-column is requested, we need another geometry column
     1241            if(Options->enable_centroid && strcmp(tables[i].type, "POINT") != 0)
     1242            {
     1243                pgsql_exec(sql_conn, PGRES_TUPLES_OK, "SELECT AddGeometryColumn('%s', 'centroid', %d, 'POINT', 2 );\n",
     1244                            tables[i].name, SRID );
     1245               
     1246                pgsql_exec(sql_conn, PGRES_COMMAND_OK, "ALTER TABLE %s ALTER COLUMN way SET NOT NULL;\n", tables[i].name);
     1247               
     1248                tables[i].hasCentroidColumn = 1;
     1249            }
     1250           
    12131251            pgsql_exec(sql_conn, PGRES_COMMAND_OK, "ALTER TABLE %s ALTER COLUMN way SET NOT NULL;\n", tables[i].name);
    12141252            /* slim mode needs this to be able to apply diffs */
    12151253            if( Options->slim )
    … …  
    12771315        if (Options->enable_hstore) strcat(sql,",tags");
    12781316
    12791317        tables[i].columns = strdup(sql);
    1280         pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN", tables[i].name, tables[i].columns);
     1318        if( tables[i].hasCentroidColumn == 1 )
     1319            pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way,centroid) FROM STDIN", tables[i].name, tables[i].columns);
     1320        else
     1321            pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN", tables[i].name, tables[i].columns);
    12811322
    12821323        tables[i].copyMode = 1;
    12831324    }

  • output-gazetteer.c

     
    878878      xnodes[count] = NULL;
    879879      xcount[count] = 0;
    880880
    881       wkt_size = build_geometry(id, xnodes, xcount, 1, 1, 1000000);
     881      wkt_size = build_geometry(id, xnodes, xcount, 1, 1, 1000000, 0);
    882882      for (i=0;i<wkt_size;i++)
    883883      {
    884884         char *wkt = get_wkt(i);

  • parse-xml2.c

     
    347347            EndElement(name, osmdata);
    348348            break;
    349349        case XML_READER_TYPE_SIGNIFICANT_WHITESPACE:
     350        case XML_READER_TYPE_COMMENT:
    350351            /* Ignore */
    351352            break;
    352353        default:
    353             fprintf(stderr, "Unknown node type %d\n", xmlTextReaderNodeType(reader));
     354            fprintf(stderr, "Unknown xml-node type %d\n", xmlTextReaderNodeType(reader));
    354355            break;
    355356    }
    356357