Dyalog APL, 186 179 bytes

I←⌷⍨∘⊃⍨⍤0 99
a←(⊂¨⍸g) {d←(⍴⍵)⍴⌊/⍬ ⋄ d[⍺]←0 ⋄ ((~⍺⍷⍸⍵)/⍸⍵)({t←⍵ ⋄ t[⍺]←⌊/↑1+({⍵+⊂(-,⊢)(⌽,⍥⊂⊢)1 0}⊂¨⍺)I¨⊂t ⋄ t}⍣≡)d}¨ ⊂g←i='#'
x←{t[⊃1↓⍒t←∪,⍵]}¨a
m←⊃x[⍒x]
i[⍸(⊃a[⊃⍒x])∊0 m]←'X'
i

De-golfed:

I←⌷⍨∘⊃⍨⍤0 99                          ⍝ Sane indexing from https://xpqz.github.io/learnapl/
neighbors←{⍵+⊂(-,⊢)(⌽,⍥⊂⊢)1 0} ⍝ Returns the cardinal neighbors of each element
                   (     )     ⍝ As a tacit function
                        ⊢      ⍝ take the array 1 0
                    ⌽          ⍝ and a reverse copy
                     ,⍥⊂       ⍝ and join them after encapsulating each of them
                               ⍝ (0 1) (1 0)
              (   )            ⍝ And as another tacit function
                 ⊢             ⍝ take that new array
               -               ⍝ and a negated copy
                ,              ⍝ and join them together
                               ⍝ (0 -1) (-1 0) (0 1) (1 0)

explore←{t←⍵ ⋄ t[⍺]←⌊/↑1+(neighbors⊂¨⍺)I¨⊂t ⋄ t} ⍝ Takes a set of coordinates of land 
                                                ⍝ on the left and a map of distances on the right.
                                                ⍝ Marks the coordinate with the minimum of the distances of the neighbors + 1

paths←{d←(⍴⍵)⍴⌊/⍬ ⋄ d[⍺]←0 ⋄ ((~⍺⍷⍸⍵)/⍸⍵)(explore⍣≡)d} ⍝ Takes a singular starting coordinate on the left
                                                      ⍝ and the grid of all land tiles on the right
       d←(⍴⍵)⍴⌊/⍬                                     ⍝ Create a grid the same shape as the input with every element set to (essentially) positive infinity for the distances map
                   d[⍺]←0                             ⍝ Set the distance at the starting coordinate to 0
                                        (explore⍣≡)   ⍝ Explore until fixed point
                            ((~⍺⍷⍸⍵)/⍸⍵)              ⍝ Remove the starting coordinate from the list of coordinates pulled from the input grid for the left argument
                                                   d  ⍝ and the to-be updated distances map as the right argument

a←(⊂¨⍸g) paths¨ ⊂g←i='#'
                   i='#'  ⍝ Set where the input is equal to a '#' to 1 and the rest to 0
                 g←       ⍝ Save the result to 'g' (for grid)
  (⊂¨⍸g)                  ⍝ Box each coordinate of a 1 in grid
         paths¨           ⍝ Call 'paths' on each of the left argument
                ⊂         ⍝ on the entirety of the grid

x←{t[⊃1↓⍒t←∪,⍵]}¨a   ⍝ Find the largest path of all
           ∪,⍵       ⍝ Ravel and uniquify the all the distances
         t←          ⍝ and store the result in t
     ⊃1↓⍒            ⍝ Sort descending and drop the first result which will always be positive infinity
   t[         ]      ⍝ Index by that to grab the actual value
  {            }¨a   ⍝ for each of the distance maps

m←⊃x[⍒x]             ⍝ Sort x descending and get the first (highest) value which is the maximum distance between two land spots
                     ⍝ This is to help find it on the distance map

i[⍸(⊃a[⊃⍒x])∊0 m]←'X' ⍝ Mark the 0 and m spot with 'X's
   (⊃a[⊃⍒x])          ⍝ Pick the best map
            ∊0 m      ⍝ Find elements corresponding to 0 and m
  ⍸                   ⍝ find their coordinates
i[              ]←'X' ⍝ and set them to 'X'

i                     ⍝ Return the final map

Python, 198 bytes

def f(s):
 m=s.index("\n")+1;x=y=z=0
 for i in range(len(s)):
  for j,k in zip(q:=[i],d:=[0]):
   if s[j]=='#':
    if k>z:x,y,z=i,j,k
    d+=[k+1]*len(e:={j+1,j-1,j+m,j-m}-{*q});q+=e
 s[x]=s[y]='X'

Attempt This Online!

Uses #, . and X

Uiua, 98 bytes

Uses #, ., and $

m←⬚0>°⊚≡↘1
⬚0+°⊚≡↘1≡⊢≡⊏¤⍖≡⊢⊢⇌.⊟⟜≡(⊙◌⊢⇌⍢(⊂:⊙⊙m⟜:▽⊃(⊡◌≡°⊂⊙◌|⊙⊙∘)+≡⊂1⊂¯.⋯1_2↯4°⊂|/↥♭,)⟜m¤)⊙¤≡⊂0⊚.=@#.

Try it!

Try it more! (more test cases, ungolfed+explanation)

Python3, 565 bytes

A little longer than the other solutions, but very fast.

import math
E=enumerate
M=[(1,0),(0,1),(-1,0),(0,-1)]
def U(i,j,d):
 q=[(*i,[i])]
 while q:
  q=sorted(q,key=lambda x:math.dist(x[2][-1],j))
  x,y,p=q.pop(0)
  if(x,y)==j:return len(p)
  for X,Y in M:
   if(T:=(x+X,y+Y))not in p and T in d:q+=[(*T,p+[T])]
def f(b):
 d={(x,y)for x,r in E(b)for y,v in E(r)if'#'==v}
 w=[(x,y)for x,y in d if any('.'==b[x+X][y+Y]for X,Y in M+[(1,-1),(1,1),(-1,-1),(-1,1)])]
 t=[]
 for I,i in E(d):
  for j in w[I+1:]:t+=[(*i,*j,U(i,j,d))]
 x,y,X,Y,L=max(t,key=lambda x:x[-1])
 b[x][y]='X'
 b[X][Y]='X'
 return'\n'.join(map(''.join,b))

Try it online!

Python 3, 249 246 bytes

Shaved off 3 bytes, thanks DLosc.

Input and output are single strings, with '.', '@', and 'X' representing water, huts, and land, respectively.

A='@'
def f(s):
 w=s.find('\n')+1;d=u={(k,k):0for k,c in enumerate(s)if A<c}
 while u:d.update(u);u={(k,j):d[(k,i)]+1for k,i in d for j in{i+1,i+w,i-1,i-w}if A<s[j]and(k,j)not in d}
 k,j=sorted(max(d,key=d.get))
 return s[:k]+A+s[k+1:j]+A+s[j+1:]

Prior version:

Input is a single string, with '.' and '#' representing water and land, respectively. 'X' represents the huts in the output.

def f(s):
 w=s.find('\n')+1;d=u={(k,k):0for k,c in enumerate(s)if'#'==c}
 while u:d.update(u);u={(k,j):d[(k,i)]+1 for k,i in d for j in{i+1,i+w,i-1,i-w}if'#'==s[j]and(k,j)not in d}
 k,j=sorted(max(d,key=d.get))
 return s[:k]+'X'+s[k+1:j]+'X'+s[j+1:]

Explanation:

It's basically doing a breadth first search from every possible starting point at the same time. Keep a dictionary, d, of path lengths keyed by the start and end of the path, e.g., d[(k,i)] is the distance from k to i. Then iterate over the keys in the dictionary, d, and create a new dictionary, u, with paths that are 1 unit longer by moving the end point 1 unit to the N,S,E,W, e.g., u[(k, i+1)] = d[(k,i)] + 1. Don't include paths that are already in d. If u is not empty, then add the new longer paths to d and repeat. When u is empty, that means no more paths could be made. Now d contains all the possible paths and their lengths. So its just a matter of getting the key with the longest path.

Less golfed, commented version:

def f(s):
  w=s.find('\n')+1                    # width of a row, or a move N or S

  d = {}                              # dictionary of all the paths.
                                      # The key is a tuple (k,j) and the
                                      # value is the distance from k to j.
  for k,c in enumerate(s):            # Initialize. Distance from k to k is 0
    if'#'==c:                         # Only do land.
      d[(k,k)] = 0

  u = d                               # dictionary of new paths. initialize it to d
                                      # so loop is entered. first d.update is
                                      # basically a NOP

  while u:                            # while there are new paths
    d.update(u)                       # add the new paths to the dict of old paths
    u={}                              #
    for k,i in d:                     # iterate over the known paths. k is the start, i is the end
      for j in{i+1,i+w,i-1,i-w}:      # iterate over squares 1 move to the E,S,W,N from i
        if'#'==s[j]and(k,j)not in d:  # if it's still land, and the path (k,j) isn't already in d,
          u[(k,j)] = d[(k,i)]+1       # then add the new path to u

  k,j=sorted(max(d,key=d.get))        # find the longest path

  return s[:k]+'X'+s[k+1:j]+'X'+s[j+1:]  # X marks the endpoints.

C#, 387 bytes

Let's get the ball rolling...

using C=System.Console;class P{static void Main(){string D="",L;int W=0,H=0,z,n,q,h,b=0,c,a,i=0,j=0;for(;(L=C.ReadLine())!=null;H+=W=L.Length)D+=L+="\n";for(z=H;z-->0;){int[]S=new int[H],Q=new int[H*8];for(Q[h=q=0]=z;q<=h;)if((c=S[n=Q[q++]]-1)<0&D[S[n]=n]==35)for(a=4;a-->0;b=c<b?c+(i=z)*(j=n)*0:b)S[Q[++h]=new[]{1,-1,W,-W}[a]+n]=S[Q[h]]<1?c:1;}for(;++z<H;)C.Write(z==i|z==j?'X':D[z]);}}

Try it Online

Complete program, reads from STDIN, writes to STDOUT. It simply goes over each cell, and runs a BFS to compute the farthest cell, recording both if it is the farthest on record. Nothing to it really, and frustratingly little I can find to golf.

Formatted and commented code:

using C=System.Console;

class P
{
    // \n 10
    // \r 13
    // . 46
    // # 35
    // x 88

    static void Main()
    {
        string D="", // map
            L; // line of input

        int W=0, // width
            H=0, // length
            z, // outer position
            n, // next position to expand
            q, // queue position pointer
            h, // queue head pointer
            b=0, // best
            c, // distance to this cell (negative)
            a, // counter
            i=0, // hermit 1 pos
            j=0; // hermit 2 pos

        for(;(L=C.ReadLine())!=null; // read a line, while we can
                H+=W=L.Length) // record the width, and add to length
            D+=L+="\n"; // add a newline, and add the line to the map

        for(z=H;z-->0;) // for each cell
        {
            int[]S=new int[H], // 'seen' >0 -> seen, else it is the distance we have found to it
                Q=new int[H*8]; // due queue (fewer than H*4 expantions, two ints each)

            // standard BFS
            for(Q[h=q=0] // reset currect 
                =z; // expand z first
                q<=h;)
                if((c=S[n=Q[q++]]-1)<0& // record 'seen', and check we havn't been seen
                    D[S[n]=n]==35) // mark as seen, and check we are a hash #
                    // 'move'
                    for(a=4;a-->0; // for each of the 4 neighbours
                            b=c<b? // if we have beaten the best
                            c+(i=z)*(j=n)*0: // set new best, record hermit positions
                            b)
                        S[Q[++h]=new[]{1,-1,W,-W}[a]+n]= // queue it for expantion
                        S[Q[h]]<1? // not seen? (this is BFS, don't need to check c is less thatn S[l+n]
                        c: // distance
                        1; // mark as seen (means it won't actually be expanded)
        }

        // z = -1
        for(;++z<H;) // for each cell
            C.Write(z==i|z==j?'X':D[z]); // print either the original char, or X if it is a hermit's home
    }
}