Thursday, June 23, 2011

Processing geometries

As said in the post about data types, geometries are simple data types at the same level strings are. Typically this means that they are treated as a whole and there is no way to mess around with the internal components of those values.

However, in case of GGL2 there is a special construction to process the geometry coordinates individually: the "process coordinates" block.
g = LINESTRING(0 0, 10 10);
process g coordinates as coordinate {
show coordinate;
}
These construction lets us process all the coordinates of a geometry, but there's more. If a "replace" statement is used inside the code block it will replace the values of that coordinate with the values specified in that same statement. For example, this makes possible to modify internally a geometry to translate all the coordinates of a geometry with an offset:
process geom coordinates as c {
x = ST_X(c) + deltaX;
y = ST_Y(c) + deltaY;
replace x y z;
}
And the previous code conveniently packaged:
alg ST_Translate(geometry geom, double deltaX, double deltaY) returns geometry {
process geom coordinates as c {
x = ST_X(c) + deltaX;
y = ST_Y(c) + deltaY;
z = ST_Z(c) + deltaZ;
replace x y z;
}
return geom;
}
This method can be used with a single geometry:
show ST_Translate(POINT(0 0), 10, 10);
but also in a "select" statement to process all the elements in a table:
show myShape select (f| ST_Translate(f/the_geom, 10, 10));
Another interesting algorithm may be the ST_Transform one, that performs coordinate system transformations on the geometries. As we'll see in the next post, it is possible to call Java code from GGL2 so it would be relatively easy to implement it using geotools... or proj4.


Tuesday, June 14, 2011

Packaging and reusing algorithms

Now that I've shown you some little tricks with GGL2 it's time to see how we can package some of these tricks for a later use.

Consider the previous post example, the grid creation:
grid = [as geometry];
width = 100;
height = 100;
foreach x1 in 0 .. 1000 step width {
foreach y1 in 0 .. 1000 step height {
double x2 = x1 + width;
double y2 = y1 + height;
geometry polygon = POLYGON((x1 y1, x2 y1, x2 y2, x1 y2, x1 y1));
add polygon to grid;
}
}
It creates a grid that covers the (0, 0)-(1000, 1000) extent with cells 100x100 big. Having such a piece of code may be useful, but having to change manually the values of the extent or the size of the cell is not handy. The solution are "algorithms":
library ggl.grid;

import ggl.geom;

alg createGrid(double cellWidth, double cellHeight, geometry extent) returns sequenceof geometry {
grid = [ as geometry ];
foreach x in ST_MinX(extent) .. ST_MaxX(extent) step cellWidth {
foreach y in ST_MinY(extent) .. ST_MaxY(extent) step cellHeight {
double x_ = x + cellWidth;
double y_ = y + cellHeight;
geometry polygon = POLYGON((x y, x_ y, x_ y_, x y_, x y));
add polygon to grid;
}
}
return grid;
}
Algorithms are contained in libraries that are declared as seen in the first statement. As we'll use some algorithms defined in the "ggl.geom" library we need to import it with the second statement. The rest of the code is similar to the previous post, but wrapped in the algorithm declaration, which defines the "variable" parts of the code as parameters so that it is possible to call the same algorithm with different values. In this concrete example, the parameters are "cellWidth", "cellHeight" and "extent".

It is possible to "call" the algorithm by naming it and passing the concrete values between parenthesis. For example, we could create a grid to cover the (0, 0)-(10, 1) extent with cells 1x1 big:
import ggl.grid;

myGrid = createGrid(1, 1, POLYGON((0 0, 10 0, 10 1, 0 1, 0 0)));
A last note, GGL2 editor can be connected to GIS instances, such as GearScape (and soon gvSIG), making it possible to reference the layers in the GIS instance by their name. In such scenario, it's very simple to create a grid that covers the whole extent of a layer:
import ggl.grid;

grid = createGrid(100, 100, ST_Extent(myLayer/the_geom));
In a next post we'll see some polymorphic features that are very useful when packaging algorithms. We'll also take a look at the possibility to call Java code from GGL2.

Monday, June 13, 2011

Create grid

Create grid

In a previous post I showed how it was possible to create a polygon using variables inside the WKT expressions. Now we'll see how it is possible to create a grid in few lines of GGL2 code.

Consider that we want to create a grid which cell size is 100x100 and that covers the (0, 0)-(1000, 1000) extent. By using two nested loops that run from 0 to 1000 in increments of 100 it is possible to obtain the upper left corner of every cell:
foreach x1 in 0 .. 1000 step 100 {
foreach y1 in 0 .. 1000 step 100 {
show x1 + ', ' + y1;
}
}
Now in the body of the inner loop we can create a cell using pseudo WKT and add it to a resulting collection:
grid = [as geometry];
width = 100;
height = 100;
foreach x1 in 0 .. 1000 step width {
foreach y1 in 0 .. 1000 step height {
double x2 = x1 + width;
double y2 = y1 + height;
geometry polygon = POLYGON((x1 y1, x2 y1, x2 y2, x1 y2, x1 y1));
add polygon to grid;
}
}
Simple, right? The only thing that may worth explaining is the creation of an empty sequence in the first line... but that will be another post. I'll also show how to package this as an algorithm and to use the extent of an existing data source as the area to cover with the grid.