Mathematica (version 9), 166 bytes

The nice, short ConvexHullMesh function that Greg Martin used was only introduced in Mathematica version 10, so I thought I'd make an attempt without it, using my ancient Mathematica version 9. I managed to get it slightly shorter, though! It's a function that takes and returns a list of strings (with ., # and o as the symbols).

f@m_:=m(1-m[[2,2]])CrossMatrix@1;""<>#&/@("o"MorphologicalComponents[#,Method->"ConvexHull"]~MorphologicalTransform~f+"#"#/.{0->"."})&[Characters@#/.{"."->0,"#"->1}]&

Explanation:

• First, Characters@# /. {"."->0, "#"->1} turns the input into a matrix of 0s and 1s.
• "o" MorphologicalComponents[#, Method->"ConvexHull"] ~MorphologicalTransform~ f + "#" # then uses Mathematica's powerful (but extremely byte-heavy…) image processing capabilities to first fill in the island's convex hull (which is the shape you'd get if you stretched a piece of string around it), and then take its boundary (using the function f@m_:=m(1-m[[2,2]])CrossMatrix@1). We then multiply this matrix by the string "o" to get a matrix of 0s and "o"s (thanks to Mathematica's impressive adaptability about types), and add this to "#" times the matrix of the island.
• Finally, ""<>#& /@ (... /. {0->"."}) turns this matrix of "o"s, "#"s and 0s into a matrix of "o"s, "#"s and "."s, and joins each row into a string.

When we test this on example B, we get the output

{"....oo..",
 "...o##o.",
 "..o####o",
 ".o###..o",
 "o#####o.",
 "o#####o.",
 ".o##oo..",
 "..oo...."}

Mathematica, 168 bytes

(x_~n~y_:=Min[Norm[x-#]&/@y];i=#;p=i~Position~#&;t=p["#"|"."]~Select~#&;(i~Part~##="o")&@@@t[#~n~t[ConvexHullMesh[Join[s=p@"#",#+{.1,.1}&/@s]]~RegionMember~#&]==1&];i)&

Pure function taking a 2D array of characters as input and returning a 2D array of characters. An easier-to-read version:

1  (x_~n~y_ := Min[Norm[x - #] & /@ y];
2  i = #; p = i~Position~# &; 
3  t = p["#" | "."]~Select~# &;
4  (i~Part~## = "o") & @@@ 
5    t[#~n~
6      t[ConvexHullMesh[
7        Join[s = p@"#", # + {.1, .1} & /@ s]]
8      ~RegionMember~# &] == 1 &];
9  i) &

Line 1 defines a function n that produces the (smallest) distance between a point x in the plane and a set y of other points. Line 2 initializes the variable i to the input, both to resolve a currying ambiguity later and to be able to modify it to produce the eventual output; line 2 also defines a function p that returns the coordinates of all occurrences of its input in i.

On line 3, p["#" | "."] represents every single coordinate from the input map (since all its characters are either "#" or "."), so t is a function that selects only those coordinates that satisfy a yet-unspecified property. On line 4, i~Part~## = "o" is going to change a bunch of entries of i to the character "o"; those characters will be selected from the set of possible coordinates according to the stuff on lines 5-8. And line 9 just returns the answer once it's computed.

Okay, infrastructure done, now to the actual computation. ConvexHullMesh is Mathematica's built-in for computing the convex hull of a set of points (the smallest convex polygon containing those points). Morally speaking, this should "fill in" the coves and fjords of the island (which is s = p@"#"), to rule them out of our cicrumnavigation. There's a little problem with ConvexHullMesh when that set of points is all in a line (thank you, test case #2), which we resolve by appending a slightly offset version of s to itself in line 7. This output is a polygon, so lines 7-9 (t[...~RegionMember~# &]) produces a list of the points with integer coordinates in that polygon. Finally, line 5 and the end of line 9 calculate all points that are at a distance exactly 1 (hence not 0) from this set of integer points; that set becomes the circumnavigating path.

Below is the output for the large test case at the OP's link. Notice at the upper left, the unusual choices of when to go west versus southwest hint at the fact that it's shadowing an invisible line of slope -2/3 between two peninsulas (said line segment being part of the convex hull's boundary).

........................
.............o..........
...........oo#ooooooo...
..........o#.#.##...#o..
........oo.#.#.###.##o..
.......o..########.##o..
.....oo...############o.
...oo#....############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#.....###....#oooo.....
.oooooo#ooooooo.........
.......o................

Python 2, 430 bytes (not competing)

It appears to fail on specific testcases. I let it be here to help others. Also I hope to fix it.

Lengthy crude solution.
It scans array from top to bottom, exept first and last row, finding most left and most right land tiles in each row.
Places circumnavigation point near finded tile, unless tiles on current row shifted inwards more than 1 from tiles from previous row. In that case circumnavigation point is placed 1 off inwards, relatively on border tiles from previous row.
After all line are scanned, similar scan from bottom to top is performed, to fixed all offset circumnavigation points. Offsets occure when row in shorter than the one under it.
Last step - close circumnavigation at first and last rows (without land) I'm too sleepy to squeeze it more, will try tomorrow.

Try it online (takes input as list of lists of characters)

L=len;R=range
S=input()
A=[]
a=0
for s in S[1:-1]:l=s.index('#');r=L(s)+~s[::-1].index('#');a=a and[min(l-1,a[0]+1),max(r+1,a[1]-1)]or[l-1,r+1];A+=[a]
for i in R(L(A)-1,0,-1):
 l=A[i];u=A[i-1]
 if u[0]-l[0]>1:A[i-1][0]=l[0]+1
 if l[1]-u[1]>1:A[i-1][1]=l[1]-1
for i in R(L(A)):S[i+1][A[i][0]]='o';S[i+1][A[i][1]]='o'
x=A[0][0]+1
y=A[0][1]
v=A[-1][0]+1
w=A[-1][1]
S[0][x:y]='o'*(y-x)
S[-1][v:w]='o'*(w-v)
print '\n'.join(map(str,S))