Imagine you were given a map divided into blocks, with a lake marked on it. And now you're supposed to decide to which block this lake belongs.
Map with blocks and lake
For purposes of this post, let's assume we should select the block with the biggest overlap with the lake itself. In this case, both the blocks and the lake are represented as
ST_GEOMETRY types within our database. And we are accessing our database using the AMDPs.
The HANA database allows native handling of spatial data (see
here for HELP page).
ST_INTERSECTS
To start us off, let's select just the blocks that have some kind of intersection with our lake. For this end, we can use the
ST_Intersects( ) method. It is available for all
ST_Geometry types.
(Help page)
If an intersection exists, this will return 1. Otherwise, 0.
If (BLOCKLOCATION.ST_Intersection(:lake) = 1)
Then
-- this Blocklocation does have an intersection with our lake
End if;
So, we can get all the blocks that have some intersection with our lake. Now we need to sort them by their size.
ST_INTERSECTION
We can get an
ST_Geometry representation of the intersection using the
ST_Intersection( ) method.
(Help page)
It returns an empty geometry if there is no intersection for the given two geometries. Otherwise, it return a geometry that describes the intersection.
Intersection = Blocklocation.st_intersection(:lake);
ST_AREA
Since we do have the geometry representation of the intersection, it should be easy to compute it's area, right? Just use
ST_Area( ) method on the intersection?
Intersection_area = intersection.ST_Area( )
Well, not so quickly. The
ST_Area method can be used with geometries of types
ST_Polygon and
ST_MultiPolygon.
What are the possible results of the
ST_Intersection method? Well, any
ST_Geometry. Which can include
ST_Point,
ST_Line and
ST_GeometryCollection. And these will cause issues. So let's look at these in detail.
ST_DIMENSION
ST_Point and ST_Line
It can happen that the whole intersection is actually just a point/line where the geometries touch. But they do not actually overlap. So, how can we filter these out? For example, by using the
ST_Dimension( ) method.
(Help page)
If (intersection_area.st_dimension( ) = 0 or intersection_area.st_dimension( ) = 1)
Then
-- these are points and lines, not interesting for us now
End if;
ST_GEOMETRYTYPE
ST_GeometryCollection
Another possibility is, that there can be multiple intersections. Like in this picture.
In this case, the value would be of type
ST_GeometryCollection. It's dimension would then be equal to the highest dimension of its parts. So, if our collection contains a line and a polygon, the dimension would be 2.
How can we recognize that we're dealing with a collection? We can use the
ST_GeometryType( ) method.
If (intersection_area.st_geometryType( ) = 'ST_GeometryCollection')
Then
-- now we know that our intersection area is a collection
End if;
So, how do we calculate the area of a collection? One option is to iterate over all the geometries within the collection and add all their respective areas.
ST_NUMGEOMETRIES
For iteration, it would be nice to know the number of geometries in the collection. Fortunately, there is a method for this, the
ST_Numgeometries method.
Help Page
And once we know the overall number, we can iterate over the collection and retrieve the individual geometries using the
ST_GeometryN method.
Help page
for index_coll in 1..intsec.st_numgeometries( )
do
-- now we can iterate over the collection and access the Nth geometry each iteration like this
Geom = intsec.st_geometryN( :index_coll );
End for;
Putting it all together
So, how do we put this all together?
Let's return to our 4 blocks and a lake from the beginning and say we have:
- Lake -- the geometry representing the lake
- Blocks -- a table that has the geometries representing our blocks (geometry), plus their numbers (blocknumber)
-- we will use theselater for iterating
DECLARE index_block, index_coll int;
-- here we get all the blocks with some intersection
Tmp_blocks = SELECT blocknumber,
b.geometry.st_intersecion(:lake) AS intersection
FROM blocks b
WHERE b.geometry.st_intersects(:lake) = 1;
-- we get the intersections with actual area and for polygons, we calculate the area straight
Tmp_areas = SELECT blocknumber,
Intersection,
Intersection.st_dimension( ) AS dimension,
Intersection.st_geometrytype( ) AS geomtype,
CASE intsec.st_geometrytype( )
WHEN 'ST_Polygon' THEN intersection.st_area( )
WHEN 'ST_MultiPolygon' THEN inertsection.st_area( )
ELSE 0
END AS intsec_area
FROM :tmp_blocks
WHERE intersection.st_dimension( ) = 2;
-- Now let's iterate over the collections and sum the areas of their parts
FOR index_block IN 1..record_count(tmp_areas)
DO
-- check that we're dealing with a collection
IF ( :tmp_areas.geomtype[ :index_block ] = 'ST_GeometryCollection' )
THEN
DECLARE intsec st_geometry = :tmp_areas.intersection[:index_block];
-- iterate over the collection
FOR index_coll IN 1..intsec.st_numgeometries( )
DO
-- check if this geometry has an area
IF (:intsec.st_geometryN( :index_coll ).st_dimension( ) = 2)
THEN
-- add it's area to the sum
tmp_areas.intsec_area[:index_block] = :tmp_areas.intsec_area[:index_block] +
:intsec.st_geometryN(:index_coll).st_area( );
END IF;
END FOR;
END FOR;
And what do we have now in our
tmp_areas table? A list of blocks that have an intersection with our lake (that is not just a border touch) and for each we have a intsec_area value, in which is the sum of the overlaps of this block with our lake. So now we just select the one with biggest value and we're done!
Conclusion
In this post, I have provided a very brief overview on a possibility how to calculate the areas of intersections of geometries using HANA native spatial data handling capabilities. I went over the basics of usage of the relevant methods, like
ST_Intersects, ST_Intersection, ST_Area, ST_GeometryType, ST_Dimension, ST_Numgeometries and
ST_GeometryN.