Working on PostgreSQL (What are join removals)

What are join removals?

Since September 2009 PostgreSQL's query planner has supported removal of LEFT JOINs. In order to properly explain what "join removals" are, we'd best look at some examples.

Let's keep it as simple as possible, say we had 2 tables like:

CREATE TABLE product (
  product_id INT NOT NULL,
  description VARCHAR(64) NOT NULL,
  PRIMARY KEY (product_id)
);

CREATE TABLE sale (
  sale_id INT NOT NULL,
  product_id INT NOT NULL,
  quantity NUMERIC(10,3) NOT NULL,
  PRIMARY KEY (sale_id)
);

If someone were to write a query such as;

SELECT s.product_id,p.description,s.quantity
FROM sale s
LEFT JOIN product p ON s.product_id = p.product_id;

The query planner would construct some sort of plan which accessed both sale and product, and joined the two tables together in order to produce the required results.

Now let's throw a crazy idea out there... What if the query joined to the product table, but the SELECT clause contained no actual columns from the product table...

SELECT s.product_id,s.quantity
FROM sale s
LEFT JOIN product p ON s.product_id = p.product_id;

Ok, perhaps it might seem strange to do this. Who would write such a query? Why bother with the join at all?

Well, if you imagine that someone had created a view which made available columns from both tables, but the consumer of the view didn't select any columns from the product table, then would the query executor really need to still make that join?

The answer is: 

It depends on:

1 We know that with LEFT JOINs, that if the product table happened not to have a matching record for a sale record that the semantic of LEFT JOIN would cause NULL values to appear in the result set for columns from the table being LEFT JOIN. 

2. Another affect that the LEFT JOIN could have would be that, if there happened to be more than 1 record matching a given product ID, then the join would match both of those rows to the sale table's row and cause the row in the sale table to be duplicated.

Item 1 would really not matter to us, since the join would only appear to be surplus to requirements none of the product table's columns were used in either the SELECT clause, WHERE clause etc. 

Item 2 really could be a problem as removing the join in if join would cause duplicate sale rows to be produced in the results would cause query results to change, which is obviously not acceptable!

If we were able to prove that these duplicate rows were not possible, then we'd be ok to remove that join. Luckily this is quite easy to prove in this case as in the example above the product_id is the primary key of the product table, therefore there can never be any duplicates product_id values in the table. So we can be certain that he two queries below will produce the same results:

SELECT s.product_id,s.quantity
FROM sale s
LEFT JOIN product p ON s.product_id = p.product_id;

and

SELECT s.product_id,s.quantity
FROM sale s;

Any although the join was explicitly asked for, the planner can simply pluck the join out of the plan and nobody should notice any sort of difference, apart from a nice gain in performance!

The work to implement this was done by Robert Haas and committed by Tom Lane back in September 2009 in commit 488d70ab46311386801c10691196ec8d755f2283

Optimising line of sight calculations with caches

In a previous blog here http://davidrowley.blogspot.co.nz/2013/08/calculating-line-of-sight-with-postgis.html I explained about some functionality I put together to allow line of sight calculations using elevation rasters along with PostGIS2.1. Before reading on you should read the previous blog to get yourself up to speed with the problems which we’re going to try to solve here.

The function which calculated the line of sight could do it fairly quickly for average line of sight lengths. If we only had to calculate 1 or 2 short lines of sight we could probably get away with the 50 milliseconds it took for the calculations to process, though if many locations need to be calculated then performance could be a problem. In tests I saw some of my queries taking up to 3 seconds when I was performing queries to find peaks within a 150 kilometre radius of a panorama.

For this problem it’s quite practical to cache the results of the calculations. I by no means have billions of point combinations which I must store line of sight results for. In my website I have around 300 panoramas and even if every panorama had 1000 mountains in a 150 kilometre radius, I’d still only have 300,000 items to cache. Caching seems ideal in this case, our elevation rasters are not going to be out of date in the near future. It’s far more likely that a change would be made to the line of sight function to make it more accurate, but in this case it would be a simple matter of purging the cache to allow the updated function to re-calculate the results.

The first step we need to take here is to create a table which can act as the cache.

CREATE TABLE line_of_sight_cache (
  viewing_location GEOGRAPHY NOT NULL,
  target_location GEOGRAPHY NOT NULL,
  isvisible BOOLEAN NOT NULL,
  PRIMARY KEY(viewing_location, target_location)
)

Here we’re just storing 2 locations along with a Boolean variable to tell us if there is a line of sight between the 2 locations.

We then need to make a few changes to the line of sight function to allow it to write to the cache once it calculates a line of sight.

CREATE OR REPLACE FUNCTION pano_is_point_visible_from(IN startloc GEOGRAPHY, 
                                                      IN endloc GEOGRAPHY,
                                                      IN writecache BOOLEAN
) RETURNS BOOLEAN AS
$BODY$
-- This function calculates if there is a line of sight between 2 points.
-- It does this by analysing elevation data along the direct path towards
-- the destination point.
-- The comments in the function describe the process as if a person was standing
-- at the startloc (parameter) and was walking towards the endloc (parameter)
-- the walking process naturally involves steps. These steps are similar to
-- what the function does. It starts by finding the elevation at both the startloc
-- and the endloc and then calculates the pitch/angle of the line of sight.
-- The function then enters a loop, this is the walking loop, where we start
-- taking steps towards the endloc. After every step we calculate the pitch from
-- the starting location to the end location. If this pitch is ever higher or equal
-- to the pitch of the endloc then we know we cannot see the endloc from the startloc
-- as it's being obscured by our current location.
--
-- The function calculates these pitches using trig functions then it takes that
-- angle which has been calculated and subtracts the number of degrees around the
-- world that the object is away. It does this over a fixed sized sphere rather than
-- a complex spheoid. Please note that at this time atmospheric refraction is not
-- taken into account, this means that objects on the horizon could be miscalculated
-- by around half a degree! This means mountains in the distance will actually be
-- visibly slightly higher than they will be said to be by this function.

DECLARE elevationdata RASTER;
DECLARE tot_distance DOUBLE PRECISION;
DECLARE cur_distance DOUBLE PRECISION;
DECLARE bearing DOUBLE PRECISION;
DECLARE start_elevation DOUBLE PRECISION;
DECLARE end_elevation DOUBLE PRECISION;
DECLARE cur_elevation DOUBLE PRECISION;
DECLARE step_size DOUBLE PRECISION;
DECLARE curloc GEOGRAPHY;
DECLARE end_pitch DOUBLE PRECISION;
DECLARE cur_pitch DOUBLE PRECISION;
BEGIN
 -- This query builds a raster which contains all of the raster data
 -- between 2 points. We use ST_Makeline to create a line between our
 -- starting point and our destination. It is quite possible that we'll
 -- find no or just partial raster data between our 2 points. Later we
 -- do test to see if we got some and return NULL if we found no data.
 SELECT ST_UNION(rast) INTO elevationdata
 FROM demelevation e
 WHERE ST_Intersects(rast, ST_Makeline(CAST(startloc AS GEOMETRY), CAST(endloc AS GEOMETRY)));
 
 -- If we found no data at all then we can quit... At this
 -- point we have no idea if there is a line of sight, so we
 -- return NULL.
 IF elevationdata IS NULL THEN
  RAISE NOTICE 'No elevation data found.';
  RETURN NULL;
 END IF;

 -- We now set the elevation of our start point and our end point.
 -- Note that there currently is a bit of a fuzz factor here and I'm
 -- adding 2 metres to both these values. This is because at least for
 -- our start value our eyes are above the ground and not on the ground, 
 -- so we'll have slightly more things in sight. For the end elevation
 -- this is not quite the case but the 2 metres was added due to the
 -- shapes of some mountains. If for example the mountain is completely
 -- flat at the top and we're standing in a location lower than it, if
 -- the summit location is marked in the middle of that flat area then
 -- we'll not be able to see it. I did not find this ideal as in reality
 -- I could see the mountain, just maybe not quite the middle of the summit,
 -- so the 2 meters being added to the end_elevation is fuzz factor, it
 -- perhaps should be more complex and be happy enough with seeing a point
 -- within X number of meters of the summit. At the time of writing this
 -- I could not decide what that value should be, so it remains like this.
 start_elevation := ST_Value(elevationdata, CAST(startloc AS GEOMETRY)) + 2;
 end_elevation := ST_Value(elevationdata, CAST(endloc AS GEOMETRY)) + 2;
 
 -- Now calculate the bearing which we must "walk" to
 -- our far off destination. 
 bearing := ST_Azimuth(startloc, endloc);


 -- A variable to save the total distance we must walk.
 tot_distance := ST_Distance(startloc, endloc);
 
 -- Set our step size, smaller steps will mean more loops and slower to calculate.
 -- Also there is no point in going in smaller steps than the detail of the raster.
 -- This should match the raster resolution or be more than it for if performance
 -- is a problem.
 
 step_size = 30; -- metres
 
 -- We must now work out the pitch in degrees of our line of
 -- sight from our current location to the destination.
 
 -- Here we use atan which will give us a pitch, or angle on a triangle, since
 -- the earth is round we need to reduce the pitch by the number of degrees
 -- around the earth that the object is. We use a fixed radius for this which
 -- is not quite accurate but it will do for now. Variations caused by
 -- Atmospheric Refraction will likely cause much more variation than the shape
 -- of the planet.
 
 end_pitch := degrees(atan((end_elevation - start_elevation) / tot_distance)) - (tot_distance / (6370986 * pi() * 2) * 360);
 
 -- We now start a loop to walk to our destination.
 -- Note that we stop checking the distance when we're
 -- within step_size to the object, as there's no point in
 -- checking if we can see the object when we're standing
 -- on top of it.

 -- First we just need to take our first step...
 cur_distance := step_size;
 
 WHILE cur_distance <= (tot_distance - step_size) LOOP
  
  -- Now work out the location of our new step based on
  -- our starting location, the current distance we've
  -- travelled and the bearing to the destination.
  curloc := ST_Project(startloc, cur_distance, bearing);
  
  -- Now let's look at the elevation of the current location.
  cur_elevation := ST_Value(elevationdata, CAST(curloc AS GEOMETRY));
  
  --RAISE NOTICE 'start_elevation = %, end_elevation = %, cur_elevation = %, cur_distance = %, bearing = %', start_elevation, end_elevation, cur_elevation, cur_distance, degrees(bearing);
  
  -- Calculate the pitch to our current location, same as we did for
  -- the destination before the loop began.
  cur_pitch := degrees(atan((cur_elevation - start_elevation) / cur_distance)) - (cur_distance / (6370986 * pi() * 2) * 360);
  
  -- Now we simply check if the pitch to from our starting
  -- point to our current location is greater or equal to
  -- the pitch we calculated from the start point to the end
  -- point. If it is then we can't see the end location due
  -- to the current point appearing to be taller from the
  -- start point.
  IF cur_pitch >= end_pitch THEN
   RAISE NOTICE 'Cannot see over object at pitch %. Pitch to destination is %, start elevation %, object elevation = %, distination elevation = %, dist from start = %, dist from end = %',
    cur_pitch, end_pitch, start_elevation, cur_elevation, end_elevation, cur_distance, tot_distance - cur_distance;
   
   -- Write to the cache if enabled.
   IF writecache = True THEN
    -- We only write if the record does not already exist.
    INSERT INTO line_of_sight_cache (viewing_location,target_location,isvisible) 
    SELECT startloc,endloc,false
    FROM (VALUES(1)) AS dual
    WHERE NOT EXISTS(SELECT 1 FROM line_of_sight_cache WHERE viewing_location = startloc AND target_location = endloc);
   END IF;
   
   RETURN FALSE;
  END IF;
  
  -- Now we can take another step then start do the whole process again...
  -- That is, providing we've not arrived at our destination yet.
  cur_distance := cur_distance + step_size;
  
 END LOOP;

 -- Write to the cache if enabled.
 IF writecache = True THEN
  -- We only write if the record does not already exist.
  INSERT INTO line_of_sight_cache (viewing_location,target_location,isvisible) 
  SELECT startloc,endloc,true
  FROM (VALUES(1)) AS dual
  WHERE NOT EXISTS(SELECT 1 FROM line_of_sight_cache WHERE viewing_location = startloc AND target_location = endloc);
 END IF;
 
 RETURN TRUE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 1000;

All that’s really changed is, I added an extra parameter named writecache which tells the function if it’s ok to write to the cache. I then added some statements just before we return true or false to cache the result in the cache table before resulting the result back to the calling function or user.

All we need to do now is modify the function which gives the visible peaks in the 150km radius from the panoramas location to try to lookup the cache. Here is a simplified version of what I'm using:

SELECT osm.name,
       osm.description,
       osm.alt,
FROM osmdata osm
LEFT OUTER JOIN line_of_sight_cache losc ON loc = losc.viewing_location AND osm.latlong = losc.target_location
WHERE ST_DWithin(osm.latlong, loc, 150000)
  AND COALESCE(losc.isvisible,pano_is_point_visible_from(loc,osm.latlong, true)) = true;

I have a query like this inside a function. “loc” is a parameter I pass to the function marking the location of the panorama. You can see that I’m performing a LEFT OUTER JOIN to attempt a lookup on the cache, if the lookup fails then losc.isvisible will be NULL, so since we’re using that column in the call to COALESCE, that value is used when there is a matching item in the cache, or if nothing is found then we execute the line of sight calculation function.

How well does it perform now you say?

Well first I’ll purge the cache to ensure we're going to hit the function each time.

test=# TRUNCATE TABLE line_of_sight_cache;

TRUNCATE TABLE

Time: 10.033 ms

test=# SELECT name,elevation FROM pano_find_osmdata_points('01010000A0E610000043C5387F1342654038F3AB3940DC45C00000000000309C40',150000);

       name       | elevation

------------------+-----------

 Mount Cook       |      3755

 Mount Malte-Brun |      3198

(2 rows)

Time: 2023.655 ms

So you can see we get our results back in around 2 seconds.

The result of executing the function once will have cached the data into the cache table, so if we query the same query again we should see our query will execute much faster this time.

test=# SELECT name,elevation FROM pano_find_osmdata_points('01010000A0E610000043C5387F1342654038F3AB3940DC45C00000000000309C40',150000);

       name       | elevation

------------------+-----------

 Mount Cook       |      3755

 Mount Malte-Brun |      3198

(2 rows)

Time: 5.312 ms

In this case we’ve managed to increase performance by around 380 times. Though calculating more results in a single pass will see that performance gap increase.

Now since I have elevation data for all of my New Zealand panoramas I think I’d like to cache results for all of the lines of sight that I’ll need. So I’ll just execute a quick query to fill the cache for me.

test=# SELECT pano_find_osmdata_points(latlong,150000) FROM pano_locations WHERE countrycode = 'NZ';

Time: 99265.579 ms

The cache was filled in just under 100 seconds!

It’s now time for me to download elevation data for my panoramas outside of New Zealand.

Calculating Line of Sight with PostGIS 2.1

For the past few years I’ve been working a lot with 360 degree panoramic photos and more recently I’ve been putting more and more effort in to get to harder to reach places to get these photos. It was only just over 1 month ago now that I decided to do a 3 day walk which was dedicated to take a single panorama!

Lots of my recent panoramas have been in mountainous regions and, I’m not sure about you, but when I’m in areas with lots of mountains I normally find myself wondering, which mountain that is over there to the North-East. It wasn’t until I discovered the existence of a mobile phone application which helps with this that I decided to try and implement something into my panoramas.

My panoramas can all be viewed on my website at www.davidrowley.co.uk/pano.php here the panoramas are viewed interactively, so you can scroll around the panorama, look up, look down, round, zoom in and out. I had already been using PostGIS 2.0 for calculating distances to other panoramas, this allowed me to find the nearest panoramas from another panorama accurately.

Once I had seen the mobile application which named the peaks I decided that I would try to do something similar in my panoramas. The platform I was using at the time was perfect to allow me to do this. I use some software called krpano, this allows my flat 360x180 degree panoramas to become interactive. It also allows many other things to be done, you can have hotspots to other panoramas and even include video in your panorama to bring things to life!

KRPano is programmed by making reference to XML files which contain setup information and code to tell the panoramic viewer what to do. For me I write these XML files dynamically using PHP and I use PHP to read details about the panorama directly from my PostgreSQL 9.2 database. I was really already half way to displaying the names of the mountain peaks in my viewer. It was just a simple matter of gathering up some peak data and writing some code to layer them onto the pano. Remember I didn’t want to have to name my peaks manually… I wanted to have this done automatically.

The first step was gathering some peak information from open street maps. I now had a fairly large list of mountains loaded into a table in my database. I then added a Geography column to the table and set a trigger to update this column from the latitude, longitude and elevation columns. With about 1 day of work my panoramic viewer was showing me mountains. I had written algorithms to calculate the pitch/angle to the top of the mountain based on the elevation of the mountain and the elevation of the panorama along with the distance between the 2. The algorithm takes into account the curvature of the earth, but still will need some more work to correct for atmospheric refraction.

With the open street map layer enabled a typical scene from my panoramic viewer now looked like this

To get this far I created a single function that sits in the database and when passed a location and radius it returns a list of mountain peaks which are within that radius, along with the bearing and pitch to that peak.

When I clicked refresh on my web browser I was really happy to see the scene above. I had learned a lot about PostGIS on this day and was fairly happy with what the outcome was. It was just a matter of hovering over the blue triangle and it would tell you the name and elevation of the peak behind it!

Now, if you look more closely at the peaks in the above image then you’ll see that some of the triangles are grouped quite closely together. Let me zoom in and give you a closer look:

Annoyingly here the mountain named “Tableland” is slightly lower and further away than the slightly higher “Gordon’s Pyramid”. So really listing Tableland was wrong here as it could not actually be seen!

This was not the only problem. I was also hitting problems were certain panorama the mountains did not line up despite me being totally accurate with their locations and elevations.

In the above image you can see that Mount Cook is much lower than it’s meant to be. The problem here lies with the panorama not being perfectly level… Something which is very hard to do in mountainous areas even after you make sure your tripod is as level as you can get it.

So it seems there’s going to be a fair bit more work to do to get this working as I want it.

Calculating Line of Sight

The first step was to gather up some elevation data. Thanks to NASA and some hard work from a dedicated man in Scotland most of the world’s elevation data is available for free here http://www.viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm

I started by downloading all the elevation data for New Zealand, which when I loaded it into the 3dem software which can be found on the viewfinderpanorama.org website I got this familiar shape.

The next step was to figure out some way to load these files into my PostgreSQL database. Luckily this is very simple as PostGIS comes with a tool to turn these elevation files into an SQL file which can be imported with psql.

I put all my .hgt into 1 directory and ran the following command:

"C:\Program Files\PostgreSQL\9.2\bin\raster2pgsql.exe" -s 4326 -I -C -M "HGT\*.hgt" -F -t 75x75 public.demelevation > elevation.sql

The above command creates a file called elevation.sql which contains all the CREATE TABLE command and imports all the data then properly constraints the tables to help PostgreSQL’s query planner do its job efficiently.

Before working with this elevation data I had no idea what a raster was. So if you’re in the same boat, then it is basically a raster when you’re talking about spatial data is a tile or array of data relating to a patch of the surface of the planet. A raster contains its anchor point and contains information about the resolution of the raster, and is generally very much like a bitmap. The remainder if the data is just an array of values. In the case of the raster’s that we are working with, this array of values is elevation data in metres.

Once I loaded the above file I could find the elevation for any point in New Zealand with the following command:

SELECT ST_VALUE(e.rast, ST_SetSRID(St_MakePoint(170.1420703, -43.594937), 4326))
FROM demelevation e
WHERE ST_Intersects(e.rast, ST_SetSRID(St_MakePoint(170.1420703, -43.594937), 4326));

What’s happening here is, we first need to find the raster tile for the location that we’re looking in, so ST_Intersects in the WHERE clause will gather all rasters with data about this location. With the resulting rasters, in this case, since we have no overlap then we only should have 1. We then ask that raster for its value for the same location. When I run this command I get the following:

test=# SELECT ST_VALUE(e.rast, ST_SetSRID(St_MakePoint(170.1420703, -43.594937), 4326))

test-# FROM demelevation e

test-# WHERE ST_Intersects(e.rast, ST_SetSRID(St_MakePoint(170.1420703, -43.594937), 4326));

 st_value

----------

     3712

(1 row)

Time: 1.499 ms

This is the approximate elevation of Mount Cook, the highest mountain in New Zealand. Also you might notice that the database returned this result to me quite quickly at around 1.5 milliseconds. This is amazingly fast and it's made possible by putting a gist index on the rast column.

So now it’s time to build a function which will calculate if there is a line of sight between 2 points.

The first step we’d need to take here is to get a raster which contains all of the rasters between 2 points. Since we can merge rasters with ST_Union we can run a query like the following.

SELECT ST_Union(rast)
FROM demelevation
WHERE ST_Intersects(rast, ST_SetSRID(ST_MakeLine(St_MakePoint(171, -43), St_MakePoint(171, -44)), 4326));

This will return a single column, single row result containing a raster with all tiles merged for a line between E171 degrees, S43 degrees and E171 degrees, S44 degrees. We now need to save that raster for processing so we can determine if there’s something along that line which will block our line of sight…

I came up with the following function which takes 2 geography parameters and returns true or false depending on if there is a line of sight between the 2 locations. Null is returns if there is no raster data and currently it's a bit undefined if there is partial raster data. Please note that the function currently does not take into account atmospheric refrations, which can have quite a big effect on objects near the horizon and may very well report things not to have a line of sight when they actually do. Currently I'm not quite sure how these should be calculated as most of the methods I've seen calculated for celestial bodies rather than objects on the face of the earth. I would imagine this also must take into account the distance of the object. If you have any ideas how I could implement this then please leave a comment.

CREATE OR REPLACE FUNCTION pano_is_point_visible_from(IN startloc GEOGRAPHY,
                                                      IN endloc GEOGRAPHY
) RETURNS BOOLEAN AS
$BODY$
-- This function calculates if there is a line of sight between 2 points.
-- It does this by analysing elevation data along the direct path towards
-- the destination point.
-- The comments in the function describe the process as if a person was standing
-- at the startloc (parameter) and was walking towards the endloc (parameter)
-- the walking process naturally involves steps. These steps are similar to
-- what the function does. It starts by finding the elevation at both the startloc
-- and the endloc and then calculates the pitch/angle of the line of sight.
-- The function then enters a loop, this is the walking loop, where we start
-- taking steps towards the endloc. After every step we calculate the pitch from
-- the starting location to the end location. If this pitch is ever higher or equal
-- to the pitch of the endloc then we know we cannot see the endloc from the startloc
-- as it's being obscured by our current location.
--
-- The function calculates these pitches using trig functions then it takes that
-- angle which has been calculated and subtracts the number of degrees around the
-- world that the object is away. It does this over a fixed sized sphere rather than
-- a complex spheoid. Please note that at this time atmospheric refraction is not
-- taken into account, this means that objects on the horizon could be miscalculated
-- by around half a degree! This means mountains in the distance will actually be
-- visibly slightly higher than they will be said to be by this function.

DECLARE elevationdata RASTER;
DECLARE tot_distance DOUBLE PRECISION;
DECLARE cur_distance DOUBLE PRECISION;
DECLARE bearing DOUBLE PRECISION;
DECLARE start_elevation DOUBLE PRECISION;
DECLARE end_elevation DOUBLE PRECISION;
DECLARE cur_elevation DOUBLE PRECISION;
DECLARE step_size DOUBLE PRECISION;
DECLARE curloc GEOGRAPHY;
DECLARE end_pitch DOUBLE PRECISION;
DECLARE cur_pitch DOUBLE PRECISION;

BEGIN
      -- This query builds a raster which contains all of the raster data
      -- between 2 points. We use ST_Makeline to create a line between our
      -- starting point and our destination. It is quite possible that we'll
      -- find no or just partial raster data between our 2 points. Later we
      -- do test to see if we got some and return NULL if we found no data.

      SELECT ST_UNION(rast) INTO elevationdata
      FROM demelevation e
      WHERE ST_Intersects(rast, ST_Makeline(CAST(startloc AS GEOMETRY), CAST(endloc AS GEOMETRY)));

      -- If we found no data at all then we can quit... At this
      -- point we have no idea if there is a line of sight, so we
      -- return NULL.

      IF elevationdata IS NULL THEN
            RAISE NOTICE 'No elevation data found.';
            RETURN NULL;
      END IF;

      -- We now set the elevation of our start point and our end point.
      -- Note that there currently is a bit of a fuzz factor here and I'm
      -- adding 2 metres to both these values. This is because at least for
      -- our start value our eyes are above the ground and not on the ground,
      -- so we'll have slightly more things in sight. For the end elevation
      -- this is not quite the case but the 2 metres was added due to the
      -- shapes of some mountains. If for example the mountain is completely
      -- flat at the top and we're standing in a location lower than it, if
      -- the summit location is marked in the middle of that flat area then
      -- we'll not be able to see it. I did not find this ideal as in reality
      -- I could see the mountain, just maybe not quite the middle of the summit,
      -- so the 2 meters being added to the end_elevation is fuzz factor, it
      -- perhaps should be more complex and be happy enough with seeing a point
      -- within X number of meters of the summit. At the time of writing this
      -- I could not decide what that value should be, so it remains like this.

      start_elevation := ST_Value(elevationdata, CAST(startloc AS GEOMETRY)) + 2;
      end_elevation := ST_Value(elevationdata, CAST(endloc AS GEOMETRY)) + 2;

      -- Now calculate the bearing which we must "walk" to
      -- our far off destination.
      bearing := ST_Azimuth(startloc, endloc);

      -- A variable to save the total distance we must walk.
      tot_distance := ST_Distance(startloc, endloc);

      -- Set our step size, smaller steps will mean more loops and slower to calculate.
      -- Also there is no point in going in smaller steps than the detail of the raster.
      -- This should match the raster resolution or be more than it for if performance
      -- is a problem.
      step_size = 30; -- metres

      -- We must now work out the pitch in degrees of our line of
      -- sight from our current location to the destination.
      -- Here we use atan which will give us a pitch, or angle on a triangle, since
      -- the earth is round we need to reduce the pitch by the number of degrees
      -- around the earth that the object is. We use a fixed radius for this which
      -- is not quite accurate but it will do for now. Variations caused by
      -- Atmospheric Refraction will likely cause much more variation than the shape
      -- of the planet.

      end_pitch := degrees(atan((end_elevation - start_elevation) / tot_distance)) - (tot_distance / (6370986 * pi() * 2) * 360);
 
      -- We now start a loop to walk to our destination.
      -- Note that we stop checking the distance when we're
      -- within step_size to the object, as there's no point in
      -- checking if we can see the object when we're standing
      -- on top of it.

      -- First we just need to take our first step...

      cur_distance := step_size;

      WHILE cur_distance <= (tot_distance - step_size) LOOP

            -- Now work out the location of our new step based on
            -- our starting location, the current distance we've
            -- travelled and the bearing to the destination.

            curloc := ST_Project(startloc, cur_distance, bearing);

            -- Now let's look at the elevation of the current location.
            cur_elevation := ST_Value(elevationdata, CAST(curloc AS GEOMETRY));

            --RAISE NOTICE 'start_elevation = %, end_elevation = %, cur_elevation = %, cur_distance = %, bearing = %',                   start_elevation, end_elevation, cur_elevation, cur_distance, degrees(bearing);

            -- Calculate the pitch to our current location, same as we did for
            -- the destination before the loop began.
            cur_pitch := degrees(atan((cur_elevation - start_elevation) / cur_distance)) - (cur_distance / (6370986 * pi() * 2) * 360);

            -- Now we simply check if the pitch to from our starting
            -- point to our current location is greater or equal to
            -- the pitch we calculated from the start point to the end
            -- point. If it is then we can't see the end location due
            -- to the current point appearing to be taller from the
            -- start point.

            IF cur_pitch >= end_pitch THEN
                  RAISE NOTICE 'Cannot see over object at pitch %. Pitch to destination is %, start elevation %, object elevation = %, distination elevation = %, dist from start = %, dist from end = %',
                        cur_pitch, end_pitch, start_elevation, cur_elevation, end_elevation, cur_distance, tot_distance - cur_distance;
                  RETURN FALSE;
            END IF;

            -- Now we can take another step then start do the whole process again...
            -- That is, providing we've not arrived at our destination yet.
            cur_distance := cur_distance + step_size;
      END LOOP;
      RETURN TRUE;
END;
$BODY$
  LANGUAGE plpgsql STABLE
  COST 1000;

We can now check if we have a line of sight by using the function like this:

SELECT pano_is_point_visible_from(
         ST_SetSrid(St_MakePoint(170.06488,-43.72071),4326), 
         ST_SetSrid(St_MakePoint(170.1420703,-43.594937),4326)
       );

Above I asked the function if it is possible to see the summit of Mount Cook from near mueller hut.
The function tells me that it is.

Now if I'm standing at the same spot, I know that in reality that I can't see Mount Tasman... So let's check if we're getting the correct results.

test=# select pano_is_point_visible_from(ST_SETSRID(St_MakePoint(170.06488,-43.72071),4326), ST_SetSrid(St_MakePoint(170.1572352,-43.5657091),4326));
NOTICE:  Cannot see over object at pitch 4.94936475078937. Pitch to destination is 4.83525382423895, start elevation 1815, object elevation = 3002, distination e
levation = 3458, dist from start = 13380, dist from end = 5384.461198757
 pano_is_point_visible_from
----------------------------
 f
(1 row)


Time: 64.390 ms

If we read the notice raised by the function we can see that we've hit land 13.38km from our starting point and the land we hit was at 3002 metres above sea level.

You can also see that the result was calculated fairly quickly at 64 milliseconds... Though if we need to calculate many points this is likely to become a problem. In the next blog I'll talk about how I optimised the line of slight calculations to cache results in a table. This is practical for my workload and improves performance around 300-400 times. http://davidrowley.blogspot.co.nz/2013/08/optimising-line-of-sight-calculations.html

Shetland Isles

So who knew that there was a small Island right in between Orkney and Shetland?

It wasn't until I was on the ferry going from Orkney to Shetland that I was looking at the map on my GPS, when I discovered there was a small island that sits pretty much right in the middle of the mainland islands of Orkney and Shetland. I had nothing that told me anything about it, just a little dot on my map with the writing "Fair Isle" next to it. I didn't have any other information and had no means of getting any, but I was able to work out that the Ferry would be sailing passed it about 4am.
So I set my alarm.

 

The sun had not quite risen yet, but you could see a little bit of light coming from the North East horizon. I didn't know if the island would be flat or have cliffs or even if anyone lived there, but when it appeared I was quite impressed with the big sea cliffs. It really set me off in a good mood for the rest of that day.

A few hours later the ferry arrived in Lerwick and I was on my way North to the most Northerly point in the UK, which is close to an area that I sailed passed a short time before while on the way back from Iceland.

I had 2 more ferry crossings before I could get to this most Northerly point. 

My mood for the day was about as laid back as you can get. I arrived at the ferry port at the North end of the mainland of Shetland to take the ferry to and island called Yell. It was about then I discover I had no cash left, so I headed south again in search of a cash point. Back at the ferry port now with some cash, some bacon and some eggs for breakfast. I relaxed in the sun and cooked and replied to a back log of emails with the free WIFI at the ferry port. I probably missed 6 or so ferries in the process.

Waiting For The Ferry To Yell in scotland

It was early afternoon before took the short 15 minute ferry ride over to Yell. I rode to the North end of the island and took another ferry over to another Island, Unst, and rode to the North side of it.

It's times like these that I envy car drivers a bit. When it comes to leaving my bike to go walking, it is never a simple operation. My camera gear always comes with me and the empty space it leaves is taken up by my motorbike clothes. It sounds simple enough to say here, but the whole operation of moving this gear over and getting changed never seems to happen as quickly as I think it should.

I was walking towards the most Northerly Lighthouse in the UK. I had no way to actually get to the lighthouse as it was built on top of a rock which is several hundred metres North of the most Northerly part of Unst. The path I followed took me around the west side of the peninsula. Puffins! Wow, I saw hundreds of Puffins. They also don't seem to be too bothered by human presense as you can get quite close without them getting scared. Of course you should always respect them and give them some room to do their nesting chores.

This place is a haven for sea birds, there are thousands of them! Some large rocks off shore are completely covered in nesting birds. Puffins have got to be the most entertaining of them to watch.

Puffins near Muckle Flugga in scotland

The long way to Shetland

It was only a little over a month ago that I was on my way back from Iceland on the ferry. I had been travelling on my motorbike around Iceland and it was now time to start the slow journey home to Scotland. It takes 2 days for this Ferry to get from Iceland to Denmark! 2 days! Not that the ferry is slow or anything, it is just a long way. I found the most depressing part of the journey, for me was that we sailed right passed Shetland, a group of islands just off the North coast of Scotland. Now... when I say sailed passed, actually the islands are in the way and the ferry has to make a small detour from a direct route to Denmark to avoid crashing into them. I woke up in the morning the day after the ferry left the Faroe Islands only to see that we were sailing right passed Muckle Flugga Lighthouse, at the most Northerly point in the UK. I picked up a mobile phone signal from Shetland and phoned my brother, for free for the first time in months. The whole experience of watching your destination country disappear over the horizon made me a little uneasy, if the ferry could have just let me off I could have been home the next day instead of in 4 days. As the islands finally disappeared and my phone signal dropped out I began to think about Shetland... despite having been born in Scotland and having lived there most of my life... the furthest North I had been was Dunnet Head, which is the most Northerly point of the mainland of the UK.... I've never been to Shetland before!
 
I really wanted to get to Shetland to see what was there, but also find out if there was some way I could, one day, get on a boat there and get myself to the Faroe Islands so that I could get the ship back to Iceland and see some of the places that I missed.

Finally I arrived home 4 days later having ridden 1000km (600 miles) from the Ferry terminal in Denmark to the Ferry port near Amsterdam in Holland and taking the overnight Ferry to Newcastle in the North of England then riding home (another 250 km (150 miles)).

Once home, I found I didn't really have much to do there... I had a day to be photographer at a friend's wedding, then another few days to finish working on all the photos. I was bored by that time and decided that it is still summer time... I shouldn't waste too much more of it at home. I decided it was time to get back on the bike again and visit that little group of Islands that I had sailed passed just a week or so before.


My home in the South West of Scotland is about as far away from Shetland as you'll get and still be in Scotland. I really didn't want to just ride passed everything in between, so I rode up the West coast of Scotland taking a week before I got on the Ferry to the Orkney Isles.

As always I didn't really have much of a plan. I knew very little about Orkney and Shetland. The whole point of going there was to find out what was there.

About a week after leaving home, I finally reached a little town called Scrabster in the North of the mainland of Scotland. Going by my map there was a ferry link from here to Orkney, then another one from Orkney to Shetland. Visiting the ferry office I found out there was a ferry in a few hours and I managed to get myself and my motorbike booked onto that ferry despite it being one of Orkney's busiest weekends. While on the Ferry I learnt that some of the UK's biggest sea cliffs are here and that the tallest sea stack.. "The old man of Hoy" is here too.... I just had to go and see it!

I spent that night in Stromness in Orkney and made some plans to get to a smaller Island called Hoy. This is Hoy, home to the "Old man of hoy".

At 11am the next day I was on Hoy with my Motorbike looking up towards the hills. Hoy seems to be the only island out of all of islands that make up Orkney that has any sort of hills. I wanted to climb to the top of Ward Hill, which is the highest hill on Orkney, but it was pointless as the cloud was down at about 100 metres and even the old man of Hoy would have been hidden. Today was just not a good day. I made an attempt to kill some time in order to give the weather a chance to improve.

I set off to explore the island, making it to the far west side, then the far east side. I felt like a fly in a jam jar. The island felt too small and I couldn't get in the right frame of mind to want to explore. Later in the day rode over a causeway onto another smaller island (Longhope) and found myself at a lighthouse. By this time I was not feeling great, I had picked up a cold or some sort of bug while at my friends wedding and it was still haunting me. I felt like I had no energy and to avoid falling asleep while riding my bike I decided to park and rest for a while. I woke up several hours later still on the seat but hunched over the petrol tank using my hands with the gloves still on as a pillow between the handle bars and the petrol tank. I was still feeling tired, but just didn't have the energy to put up my tent. I was also worried that if I fell asleep on the bike again, I'd end up falling over and the bike landing on top of me in my sleep. This time instead of staying on the bike, I lay down on the ground and used the front wheel as a pillow and slept for another few hours.

It was now 8pm, still light. I had been sleeping for 4 hours total. I quickly put up my tent and crawled inside.

I woke the next day at 11am having slept for about 21 hours. I felt better but still not 100%. The good news was that the cloud had lifted in the night and the sun was out. After packing up I made my way over to Rackwick, where the walking path starts for the old man of hoy.

One hour later I was sitting at the top of a 100 meter cliff watching some climbers climbing to the top of the old man. I watched them for a few hours and in the meantime I managed to work out the path they had taken down the cliff to get to the base. Once the climbers made it to the top and were on their way back down I made my way down the cliff to set up my camera on the tripod and get a panoramic photo of the scene. By this time the sun had come around to the west a bit more and lowered in the sky. The sun was lighting up the sea cliffs and the sandstone walls of the old man.


The Old Man Of Hoy in scotland

 The Old Man Of Hoy in Scotland

The following morning I made my way back to the mainland of Orkney and explored some more. I learnt about the Churchill Barriers and how the British naval fleet had used Orkney as a base during World War 2. I visited a small chapel which was built by some Italian prisoners of war.


Italian Chapel in scotland

Evening came and it got to that time of day that I felt like I should find a place to sleep for the night. I glanced at the Ferry timetable that I had picked up while waiting for the Ferry to Orkney and discovered that I could get the ferry to Shetland either tonight or I'd have to wait for 2 more days. Not really knowing what the rest of Orkney had in store for me and looking at the calendar and seeing that time was not really on my side, I decided to try to get myself on the Shetland Ferry which was due to arrive in Kirkwall just before midnight.

On arriving at the Ferry terminal I discovered that due to the busy weekend of events that were on in Orkney the ferry was fully booked. They said I could go on standby and if there was any space then I could travel.... I agreed to risk it and I waited around, the Ferry arrived about 11pm. It was very dark by this time and I wondered where I could sleep if I didn't get on the ferry. I watched as the ferry was loaded, having no idea if I would be getting on it or not. I had always thought that there was always space for motorbikes, they take up so little room! But these thoughts had been crushed when I was turned away from the "full" ferry to Skye only 1 week before.

I had no idea if I'd be on Shetland the following day and no idea what to do if I had to stay on Orkney...