alg extrude(FeatureCollection buildings) returns FeatureCollection {
double height = 20;
geometry prev;
FeatureCollection ret;
foreach building in buildings {
process building/the_geom coordinates as g do {
if (prev is not null) {
x0 = ST_X(prev);
y0 = ST_Y(prev);
x1 = ST_X(g);
y1 = ST_Y(g);
face = POLYGON((x0 y0 0, x0 y0 height, x1 y1 height, x1 y1 0, x0 y0 0));
add { id = autonumeric(), the_geom = face } to ret;
}
prev = g;
}
prev = null;
}
return ret;
}
As we can see, our algorithm will receive a feature collection with 2D buildings and will return another feature collection. The returning collection will have one 3D building face for each edge in the original collection. This is, if we provide a square building, we will receive four 3D walls.
First, we set the height for the buildings (constant for all buildings) and declare the returning collection ret. Also, we declare a prev geometry variable where we will store a coordinate we need for the calculations. We will see it later. double height = 20;
geometry prev;
FeatureCollection ret;
Next, we iterate over each building in the collection and, for each building, we iterate over each coordinate
foreach building in buildings {
process building/the_geom coordinates as g do {
if (prev is not null) {
x0 = ST_X(prev);
y0 = ST_Y(prev);
x1 = ST_X(g);
y1 = ST_Y(g);
face = POLYGON((x0 y0 0, x0 y0 height, x1 y1 height, x1 y1 0, x0 y0 0));
add { id = autonumeric(), the_geom = face } to ret;
}
prev = g;
}
prev = null;
}
With this last process statement, what we do is to process each coordinate in building/the_geom. Then, inside the process expression, let's ignore the if statement for a moment and focus on this:
prev = g;
Each time we process a coordinate inside each building, we store that coordinate in the prev variable. So, if we use the current coordinate (g) and the previous coordinate (prev), we can have each edge in the original collection. All we need to do is to extrude that edge and add it to the returning collection.
Note also that once we have finished processing the coordinates of a building, we set prev to null again so the is null check is true for the first coordinate in the next building.
Note also that once we have finished processing the coordinates of a building, we set prev to null again so the is null check is true for the first coordinate in the next building.
And finally, let's take a look at the if statement, where the magic happens:
if (prev is not null) {
x0 = ST_X(prev);
y0 = ST_Y(prev);
x1 = ST_X(g);
y1 = ST_Y(g);
face = POLYGON((x0 y0 0, x0 y0 height, x1 y1 height, x1 y1 0, x0 y0 0));
add { id = autonumeric(), the_geom = face } to ret;
}
First, we check that prev is not null. Since we cannot extrude a wall from a single point, it will be null for the first coordinate on each building. Then, we get the x and y coordinates for the current and previous edge points and create a geometry using WKT. Finally, we add a new feature with a unique id (the autonumeric algorithm returns a different integer for each call) and the extruded geometry to the returning collection.
And now, here is an example of how to use this algorithm with some data. We used the Vancouver building footprints from [1]:
And now, here is an example of how to use this algorithm with some data. We used the Vancouver building footprints from [1]:
import ggl.shp;
import ggl.util;
import org.examples.extrude;
read SHP 'buildings/buildings.shp' to shp;
buildings = shp select(the_geom = shp/the_geom, id = autonumeric());
faces = extrude(buildings);
write SHP faces to '/tmp/faces.shp';And if we visualize the result in gvSIG:
It was not that difficult, ¿right? Now you can try to use a layer with different heights for each building, or even create buildings with ceilings ;-) Also, if you have questions, suggestions or any other applications for this algorithm, we will be glad to hear on the project mailing lists.
-----
[1] http://data.vancouver.ca/datacatalogue/index.htm
