C# 7 - 414 369 bytes

Edit: switched to 1D looping, computing i and j on the fly

using C=System.Console;class P{static void Main(){string D="",L;int W=0,H=0,i,j,z,k,q;for(;(L=C.ReadLine())!=null;H+=W=L.Length)D+=L+="\n";int[]P()=>new[]{i=z%W,j=z/W,-i,-j,i+j,j-i,-i-j,i-j};int[]B=null;for(z=H;z-->0;)if(D[z]%7<1)for(B=B??P(),k=8;k-->0;)if(B[k]<(q=P()[k]+1))B[k]=q;for (;++z<H;C.Write(q>8?'o':D[z]))for(q=k=8;k-->0;)q*=B[k]<P()[k]?0:B[k]==P()[k]?2:1;}}

Complete program, takes input to standard in, prints it to standard out, uses #, ., and o. For each cell, it computes a 'profile' (which is the distance over 8 directions), and records a maximum of each of these. It then writes out the whole map again, and replaces any cell which is both on a boundary and not outside of any with a 'o'. The commented code below explains how it all works.

As per my answer to Save the Geese from Extinction, this produces the smallest octagon (valid circumnavigation with largest area) which bounds the island.

Note: that for once in my life I'm using something from the current decade, and this code requires C# 7 to compile. If you do not have C# 7, there is one line that will need to be replaced, which is clearly marked in the code.

Example usage and output:

type t7.txt | IslandGolf1.exe

.........ooooooooooo....
........o....#......o...
.......o...#.#.##...#o..
......o....#.#.###.##.o.
.....o....########.##..o
....o.....############.o
...o.#....############.o
..o#.###.##############o
.o##.##################o
o.####################.o
o..##################..o
o.##################...o
o...################...o
o###################...o
o#####################.o
o.##################..o.
o####################o..
o#...##############.o...
o##...#############o....
o#.....###....#...o.....
.o.....#.........o......
..ooooooooooooooo.......

Formatted and commented code:

using C=System.Console;

class P
{
    static void Main()
    {
        // \n 10
        // # 35
        // . 46
        // o 111

        string D="", // the whole map
            L; // initally each line of the map, later each line of output
    
        int W=0, // width
            H=0, // length (width * height)
            i, // x-index
            j, // y-index
            z, // position in map (decomposed into i and j by and for P)
            k, // bound index
            q; // bound distance, and later cell condition (0 -> outside, 8 -> inside, >8 -> on boudary)
    
        for(;(L=C.ReadLine())!=null; // read a line, while we can
                H+=W=L.Length) // record the width, and increment height
            D+=L+="\n"; // add a \n to the line (the rest of the code treats this as a . cell), and add the line to the map
    
        // create profile for point
        // converts 1d to 2d very cheaply
        // each entry describes the distance in one of the 8 directions: we want to maximise these to find the 'outer bounds'
        // these 8 bounds describe 8 lines, together an octogen
        int[]P()=>new[]{i=z%W,j=z/W,-i,-j,i+j,j-i,-i-j,i-j}; // new C#7 local method syntax (if you don't have C#7, you can test this code with the line below instead)
        //z=0;System.Func<int[]>P=()=>new[]{i=z%W,j=z/W,-i,-j,i+j,j-i,-i-j,i-j}; // old lambda syntax (must pre-assign z to make static checker happy)
    
        int[]B=null; // our current bounds, initially null (must only call P() when on a #)
    
        for(z=H;z-->0;) // for each cell
            if(D[z]%7<1) // if this cell is #
                for(B=B??P(), // init B when we first hit a #
                    k=8;k-->0;) // for each bound
                    if(B[k]<(q=P()[k]+1))B[k]=q; // update bound if necessary (add one so that we define the bound _outside_ the hashes)
    
        // z=-1
        for (;++z<H; // for each cell
               C.Write(q>8?'o':D[z])) // print the cell (if q > 8, then we are on the bounds, otherwise, spit out whatever we were before)
            // check we are not 'outside' any of the bounds, and that we are 'on' atleast one of them
            for(q=k=8;k-->0;) // for each bound
                q*=B[k]<P()[k]?0: // outside bound (q=0)
                    B[k]==P()[k]?2: // on bound (if q != 0, then q becomes > 8)
                    1; // inside (preserve q)
    }
}