Hi
A while ago we talked about adding support for a centroid column to
osm2pgsql [1]. Attached you'll find a patch that adds this
functionality. I did not test it extensively, especially not with
invalid geometries, but for most cases it should work.
Unfortunately the server I tested on has only geos v3.1 and to use
interpolate on linestrings at least geos 3.2 ist required, so the
centroid point does not sit in the middle of the linesting but ..
somewhere else on the line.
Please, help me with testing and reviewing this patch.
Peter
[1]
<http://lists.openstreetmap.org/pipermail/dev/2010-September/020536.html>
Index: output.h
===================================================================
--- output.h (revision 23282)
+++ output.h (working copy)
@@ -29,7 +29,8 @@
const char *expire_tiles_filename; /* File name to output expired tiles
list to */
int enable_hstore; /* add an additional hstore column with objects key/value
pairs */
int enable_multi; /* Output multi-geometries intead of several simple
geometries */
- char** hstore_columns; /* list of columns that should be written into their
own hstore column */
+ int enable_centroid; /* Generate an additional geometry column that contains
the center of the geometry */
+ const char** hstore_columns; /* list of columns that should be written into
their own hstore column */
int n_hstore_columns; /* number of hstore columns */
};
Index: build_geometry.cpp
===================================================================
--- build_geometry.cpp (revision 23282)
+++ build_geometry.cpp (working copy)
@@ -21,6 +21,7 @@
*/
#include <iostream>
+#include <cstdio>
#include <cstring>
#include <cstdlib>
@@ -59,6 +60,7 @@
typedef std::auto_ptr<Geometry> geom_ptr;
static std::vector<std::string> wkts;
+static std::vector<std::string> wkt_centroids;
static std::vector<double> areas;
@@ -80,6 +82,9 @@
std::auto_ptr<LinearRing>
shell(gf.createLinearRing(coords.release()));
geom = geom_ptr(gf.createPolygon(shell.release(), new
std::vector<Geometry *>));
geom->normalize(); // Fix direction of ring
+
+ // this function seems not to be used anymore so centeroid
functionality is not implemented here
+ // see get_wkt_split
} else {
if (coords->getSize() < 2)
return NULL;
@@ -103,13 +108,13 @@
}
}
-
-size_t get_wkt_split(osmNode *nodes, int count, int polygon, double split_at) {
+size_t get_wkt_split(osmNode *nodes, int count, int polygon, double split_at,
int enable_centeroid) {
GeometryFactory gf;
std::auto_ptr<CoordinateSequence>
coords(gf.getCoordinateSequenceFactory()->create(0, 2));
double area;
WKTWriter wktw;
size_t wkt_size = 0;
+ char centroid[100];
try
{
@@ -129,6 +134,15 @@
std::string wkt = wktw.write(geom.get());
wkts.push_back(wkt);
areas.push_back(area);
+
+ // centroid of a closed way
+ if(enable_centeroid)
+ {
+ Point* p = geom->getInteriorPoint();
+ snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)",
p->getX(), p->getY());
+ wkt_centroids.push_back(centroid);
+ }
+
wkt_size++;
} else {
if (coords->getSize() < 2)
@@ -146,6 +160,23 @@
std::string wkt = wktw.write(geom.get());
wkts.push_back(wkt);
areas.push_back(0);
+
+ // centroid of a linear way
+ if(enable_centeroid)
+ {
+ // the interior point is always on an edge of the line
+ Point* p = geom->getInteriorPoint();
+
+ // the centroid is not always part of the line
+ // Point* p = geom->getCentroid();
+
+ // interpolate a point on the geometry, needs geos >=
3.2
+ // Point* p = geom->interpolate(0.5);
+
+ snprintf(centroid, sizeof(centroid), "POINT(%.15g
%.15g)", p->getX(), p->getY());
+ wkt_centroids.push_back(centroid);
+ }
+
wkt_size++;
distance=0;
segment =
std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0,
2));
@@ -181,6 +212,21 @@
return result;
}
+char * get_wkt_centroid(size_t index)
+{
+ if( wkt_centroids.size() <= index )
+ return NULL;
+
+ char *result;
+ result = (char*) std::malloc( wkt_centroids[index].length() + 1);
+
+ // At least give some idea of why we about to seg fault
+ if (!result) std::cerr << std::endl << "Unable to allocate memory: " <<
(wkt_centroids[index].length() + 1) << std::endl;
+
+ std::strcpy(result, wkt_centroids[index].c_str());
+ return result;
+}
+
double get_area(size_t index)
{
return areas[index];
@@ -190,6 +236,7 @@
{
wkts.clear();
areas.clear();
+ wkt_centroids.clear();
}
static int coords2nodes(CoordinateSequence * coords, struct osmNode ** nodes) {
@@ -290,12 +337,13 @@
return 1;
}
-size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int
make_polygon, int enable_multi, double split_at) {
+size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int
make_polygon, int enable_multi, double split_at, int enable_centeroid) {
size_t wkt_size = 0;
std::auto_ptr<std::vector<Geometry*> > lines(new std::vector<Geometry*>);
GeometryFactory gf;
geom_ptr geom;
-
+ char centroid[100];
+
try
{
for (int c=0; xnodes[c]; c++) {
@@ -358,6 +406,15 @@
std::string wkt = writer.write(geom.get());
wkts.push_back(wkt);
areas.push_back(0);
+
+ // centroid of a route-relation (linear)
+ if(enable_centeroid)
+ {
+ Point* p = geom->getInteriorPoint();
+ snprintf(centroid, sizeof(centroid), "POINT(%.15g
%.15g)", p->getX(), p->getY());
+ wkt_centroids.push_back(centroid);
+ }
+
wkt_size++;
distance=0;
segment =
std::auto_ptr<CoordinateSequence>(gf.getCoordinateSequenceFactory()->create(0,
2));
@@ -456,6 +513,16 @@
std::string text = writer.write(multipoly.get());
wkts.push_back(text);
areas.push_back(multipoly->getArea());
+
+ // centroid of a multipolygon relation (area) with more
+ // then one outer way in in multi-geometry mode
+ if(enable_centeroid)
+ {
+ Point* p = multipoly->getInteriorPoint();
+ snprintf(centroid, sizeof(centroid), "POINT(%.15g %.15g)",
p->getX(), p->getY());
+ wkt_centroids.push_back(centroid);
+ }
+
wkt_size++;
}
else
@@ -466,6 +533,15 @@
std::string text = writer.write(poly);
wkts.push_back(text);
areas.push_back(poly->getArea());
+
+ // centroid of a multipolygon relation (area)
+ if(enable_centeroid)
+ {
+ Point* p = poly->getInteriorPoint();
+ snprintf(centroid, sizeof(centroid), "POINT(%.15g
%.15g)", p->getX(), p->getY());
+ wkt_centroids.push_back(centroid);
+ }
+
wkt_size++;
delete(poly);
}
Index: osm2pgsql.c
===================================================================
--- osm2pgsql.c (revision 23282)
+++ osm2pgsql.c (working copy)
@@ -441,10 +441,11 @@
EndElement(name);
break;
case XML_READER_TYPE_SIGNIFICANT_WHITESPACE:
+ case XML_READER_TYPE_COMMENT:
/* Ignore */
break;
default:
- fprintf(stderr, "Unknown node type %d\n",
xmlTextReaderNodeType(reader));
+ fprintf(stderr, "Unknown xml-node type %d\n",
xmlTextReaderNodeType(reader));
break;
}
@@ -556,6 +557,7 @@
fprintf(stderr, " \tthat start with the specified
string, eg --hstore-column \"name:\" will\n");
fprintf(stderr, " \tproduce an extra hstore column
that contains all name:xx tags\n");
fprintf(stderr, " -G|--multi-geometry\t\tGenerate multi-geometry
features in postgresql tables.\n");
+ fprintf(stderr, " -n|--centroid\t\tGenerate an additional geometry
column that contains the center of the geometry.\n");
fprintf(stderr, " -h|--help\t\tHelp information.\n");
fprintf(stderr, " -v|--verbose\t\tVerbose output.\n");
fprintf(stderr, "\n");
@@ -651,6 +653,7 @@
int expire_tiles_zoom_min = -1;
int enable_hstore = 0;
int enable_multi = 0;
+ int enable_centroid = 0;
const char *expire_tiles_filename = "dirty_tiles";
const char *db = "gis";
const char *username=NULL;
@@ -701,6 +704,7 @@
{"hstore", 0, 0, 'k'},
{"hstore-column", 1, 0, 'z'},
{"multi-geometry", 0, 0, 'G'},
+ {"centroid", 0, 0, 'n'},
{0, 0, 0, 0}
};
@@ -744,6 +748,7 @@
hstore_columns[n_hstore_columns-1] = optarg;
break;
case 'G': enable_multi=1; break;
+ case 'n': enable_centroid=1; break;
case 'h': long_usage_bool=1; break;
case '?':
default:
@@ -814,6 +819,7 @@
options.expire_tiles_zoom_min = expire_tiles_zoom_min;
options.expire_tiles_filename = expire_tiles_filename;
options.enable_multi = enable_multi;
+ options.enable_centroid = enable_centroid;
options.enable_hstore = enable_hstore;
options.hstore_columns = hstore_columns;
options.n_hstore_columns = n_hstore_columns;
Index: build_geometry.h
===================================================================
--- build_geometry.h (revision 23282)
+++ build_geometry.h (working copy)
@@ -32,11 +32,12 @@
int parse_wkt(const char * wkt, struct osmNode *** xnodes, int ** xcount, int
* polygon);
char *get_wkt_simple(struct osmNode *, int count, int polygon);
-size_t get_wkt_split(struct osmNode *, int count, int polygon, double
split_at);
+size_t get_wkt_split(struct osmNode *, int count, int polygon, double
split_at, int enable_centeroid);
char* get_wkt(size_t index);
+char * get_wkt_centroid(size_t index);
double get_area(size_t index);
-size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int
make_polygon, int enable_multi, double split_at);
+size_t build_geometry(int osm_id, struct osmNode **xnodes, int *xcount, int
make_polygon, int enable_multi, double split_at, int enable_centeroid);
void clear_wkts();
#ifdef __cplusplus
Index: output-pgsql.c
===================================================================
--- output-pgsql.c (revision 23282)
+++ output-pgsql.c (working copy)
@@ -51,6 +51,7 @@
unsigned int buflen;
int copyMode;
char *columns;
+ int hasCentroidColumn;
} tables [] = {
{ name: "%s_point", type: "POINT" },
{ name: "%s_line", type: "LINESTRING"},
@@ -265,7 +266,11 @@
/* Return to copy mode if we dropped out */
if( !tables[table].copyMode )
{
- pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN",
tables[table].name, tables[table].columns);
+ if( tables[table].hasCentroidColumn == 1 )
+ pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way,centroid)
FROM STDIN", tables[table].name, tables[table].columns);
+ else
+ pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN",
tables[table].name, tables[table].columns);
+
tables[table].copyMode = 1;
}
/* If the combination of old and new data is too big, flush old data */
@@ -631,7 +636,7 @@
-static void write_wkts(int id, struct keyval *tags, const char *wkt, enum
table_id table)
+static void write_wkts(int id, struct keyval *tags, const char *wkt, enum
table_id table, const char *wkt_centroid)
{
static char *sql;
@@ -674,6 +679,14 @@
sprintf(sql, "SRID=%d;", SRID);
copy_to_table(table, sql);
copy_to_table(table, wkt);
+
+ if(wkt_centroid && tables[table].hasCentroidColumn == 1)
+ {
+ sprintf(sql, "\tSRID=%d;", SRID);
+ copy_to_table(table, sql);
+ copy_to_table(table, wkt_centroid);
+ }
+
copy_to_table(table, "\n");
}
@@ -800,12 +813,16 @@
else
split_at = 100 * 1000;
- wkt_size = get_wkt_split(nodes, count, polygon, split_at);
+ wkt_size = get_wkt_split(nodes, count, polygon, split_at,
Options->enable_centroid);
for (i=0;i<wkt_size;i++)
{
char *wkt = get_wkt(i);
+ char *wkt_centroid = NULL;
+ if(Options->enable_centroid)
+ wkt_centroid = get_wkt_centroid(i);
+
if (wkt && strlen(wkt)) {
/* FIXME: there should be a better way to detect polygons */
if (!strncmp(wkt, "POLYGON", strlen("POLYGON")) || !strncmp(wkt,
"MULTIPOLYGON", strlen("MULTIPOLYGON"))) {
@@ -816,15 +833,16 @@
snprintf(tmp, sizeof(tmp), "%f", area);
addItem(tags, "way_area", tmp, 0);
}
- write_wkts(id, tags, wkt, t_poly);
+ write_wkts(id, tags, wkt, t_poly, wkt_centroid);
} else {
expire_tiles_from_nodes_line(nodes, count);
- write_wkts(id, tags, wkt, t_line);
+ write_wkts(id, tags, wkt, t_line, wkt_centroid);
if (roads)
- write_wkts(id, tags, wkt, t_roads);
+ write_wkts(id, tags, wkt, t_roads, wkt_centroid);
}
}
free(wkt);
+ if(wkt_centroid) free(wkt_centroid);
}
clear_wkts();
@@ -1021,7 +1039,7 @@
else
split_at = 100 * 1000;
- wkt_size = build_geometry(id, xnodes, xcount, make_polygon,
Options->enable_multi, split_at);
+ wkt_size = build_geometry(id, xnodes, xcount, make_polygon,
Options->enable_multi, split_at, Options->enable_centroid);
if (!wkt_size) {
resetList(&tags);
@@ -1032,7 +1050,11 @@
for (i=0;i<wkt_size;i++)
{
char *wkt = get_wkt(i);
+ char *wkt_centroid = NULL;
+ if(Options->enable_centroid)
+ wkt_centroid = get_wkt_centroid(i);
+
if (wkt && strlen(wkt)) {
expire_tiles_from_wkt(wkt, -id);
/* FIXME: there should be a better way to detect polygons */
@@ -1043,11 +1065,11 @@
snprintf(tmp, sizeof(tmp), "%f", area);
addItem(&tags, "way_area", tmp, 0);
}
- write_wkts(-id, &tags, wkt, t_poly);
+ write_wkts(-id, &tags, wkt, t_poly, wkt_centroid);
} else {
- write_wkts(-id, &tags, wkt, t_line);
+ write_wkts(-id, &tags, wkt, t_line, wkt_centroid);
if (roads)
- write_wkts(-id, &tags, wkt, t_roads);
+ write_wkts(-id, &tags, wkt, t_roads, wkt_centroid);
}
}
free(wkt);
@@ -1083,11 +1105,15 @@
// If we are making a boundary then also try adding any relations which
form complete rings
// The linear variants will have already been processed above
if (make_boundary) {
- wkt_size = build_geometry(id, xnodes, xcount, 1,
Options->enable_multi, split_at);
+ wkt_size = build_geometry(id, xnodes, xcount, 1,
Options->enable_multi, split_at, Options->enable_centroid);
for (i=0;i<wkt_size;i++)
{
char *wkt = get_wkt(i);
+ char *wkt_centroid = NULL;
+ if(Options->enable_centroid)
+ wkt_centroid = get_wkt_centroid(i);
+
if (strlen(wkt)) {
expire_tiles_from_wkt(wkt, -id);
/* FIXME: there should be a better way to detect polygons */
@@ -1098,7 +1124,7 @@
snprintf(tmp, sizeof(tmp), "%f", area);
addItem(&tags, "way_area", tmp, 0);
}
- write_wkts(-id, &tags, wkt, t_poly);
+ write_wkts(-id, &tags, wkt, t_poly, wkt_centroid);
}
}
free(wkt);
@@ -1205,6 +1231,18 @@
pgsql_exec(sql_conn, PGRES_COMMAND_OK, "%s", sql);
pgsql_exec(sql_conn, PGRES_TUPLES_OK, "SELECT
AddGeometryColumn('%s', 'way', %d, '%s', 2 );\n",
tables[i].name, SRID, tables[i].type );
+
+ // if a centroid-column is requested, we need another geometry
column
+ if(Options->enable_centroid && strcmp(tables[i].type, "POINT") !=
0)
+ {
+ pgsql_exec(sql_conn, PGRES_TUPLES_OK, "SELECT
AddGeometryColumn('%s', 'centroid', %d, 'POINT', 2 );\n",
+ tables[i].name, SRID );
+
+ pgsql_exec(sql_conn, PGRES_COMMAND_OK, "ALTER TABLE %s ALTER
COLUMN way SET NOT NULL;\n", tables[i].name);
+
+ tables[i].hasCentroidColumn = 1;
+ }
+
pgsql_exec(sql_conn, PGRES_COMMAND_OK, "ALTER TABLE %s ALTER
COLUMN way SET NOT NULL;\n", tables[i].name);
/* slim mode needs this to be able to apply diffs */
if( Options->slim )
@@ -1272,8 +1310,12 @@
if (Options->enable_hstore) strcat(sql,",tags");
tables[i].columns = strdup(sql);
- pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN",
tables[i].name, tables[i].columns);
+ if( tables[i].hasCentroidColumn == 1 )
+ pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way,centroid)
FROM STDIN", tables[i].name, tables[i].columns);
+ else
+ pgsql_exec(sql_conn, PGRES_COPY_IN, "COPY %s (%s,way) FROM STDIN",
tables[i].name, tables[i].columns);
+
tables[i].copyMode = 1;
}
free(sql);
Index: output-gazetteer.c
===================================================================
--- output-gazetteer.c (revision 23282)
+++ output-gazetteer.c (working copy)
@@ -744,7 +744,7 @@
xnodes[count] = NULL;
xcount[count] = 0;
- wkt_size = build_geometry(id, xnodes, xcount, 1, 1, 1000000);
+ wkt_size = build_geometry(id, xnodes, xcount, 1, 1, 1000000, 0);
for (i=0;i<wkt_size;i++)
{
char *wkt = get_wkt(i);
_______________________________________________
dev mailing list
dev@openstreetmap.org
http://lists.openstreetmap.org/listinfo/dev
