Compute intersection size with HANA spatial data

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.

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.

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.

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.

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.

And once we know the overall number, we can iterate over the collection and retrieve the individual geometries using the ST_GeometryN method.

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!