C (gcc), 154 149 136 133 132 129 bytes

b[99];s,m;r(i){i=i>m-1u||b[i]++?0:r(i-1)+r(i+1)+r(i-s)-~r(i+s);}f(a,w,h){m=h*=s=w;for(w=0;h--;w=fmax(w,r(h)))wmemcpy(b,a,m);a=w;}

Try it online!

Recursive flood fill.

-5 -8 -9 -12 bytes thanks to @ceilingcat!

Wolfram Language (Mathematica), 94 bytes

Max[0,KeyDrop[0]@Counts@Flatten@MorphologicalComponents[#/.{}->{{0}},CornerNeighbors->False]]&

Try it online!

Scala 3, 283 bytes

283 bytes. It can be golfed more.

Golfed version. Attempt This Online!

object Main{def main(args:Array[String]):Unit={val i=scala.io.StdIn.readLine();var c=0;var G=List[Set[Int]]();for(ch<-i){var n=List(c+1);for(g<-G)if(Set(c-i.indexOf(","),c).intersect(g).nonEmpty)n=n++g;G=G:+(if("," > ch.toString)n.toSet else Set());c+=1};println(G.map(_.size).max)}}

Ungolfed version. Attempt This Online!

object Main {
  def main(args: Array[String]): Unit = {
    // Get input string (in Scala, we need to explicitly read from console)
    val inputString = scala.io.StdIn.readLine()
    
    // Initialize variables
    var currentPosition = 0
    var groups = List[Set[Int]]()
    
    // Iterate through each character in the input string
    for (character <- inputString) {
      // Check if the character is before comma in lexicographic ordering
      val isBeforeComma = "," > character.toString
      
      // Find position of first comma
      val commaPosition = inputString.indexOf(",")
      
      // Calculate position relative to comma
      val positionRelativeToComma = currentPosition - commaPosition
      
      // Create set with current position and position relative to comma
      val positionsSet = Set(positionRelativeToComma, currentPosition)
      
      // Initialize new set with current_position + 1
      var newElements = List(currentPosition + 1)
      
      // For each existing group
      for (group <- groups) {
        // Check if there's any overlap between positionsSet and the group
        if (positionsSet.intersect(group).nonEmpty) {
          // If there's overlap, add all elements from the group
          newElements = newElements ++ group
        }
      }
      
      // Create new set and multiply by isBeforeComma (true/false)
      // If isBeforeComma is false, the set will be empty
      val newGroup = if (isBeforeComma) newElements.toSet else Set[Int]()
      
      // Add the new group to our collection
      groups = groups :+ newGroup
      
      // Move to next position
      currentPosition += 1
    }
    
    // Find the maximum size of any group
    val maxGroupSize = groups.map(_.size).max
    println(maxGroupSize)
  }
}

BQN, 36 bytes^SBCS

{0⌈´1↓/⁼⥊(⍉∘⌽××⊢⌈»)⍟(4×≠⥊𝕩)𝕩×1+⊒⌾⥊𝕩}

Run online!

Replace ones in the matrix with positive integers, replace by largest integer on island by flood-filling, count each positive integer, get the largest count.

𝕩           𝕩×1+⊒⌾⥊𝕩    (⍉∘⌽××⊢⌈»)⍟(4×≠⥊𝕩) 1↓/⁼⥊            0⌈´
┌─          ┌─          ┌─                 ⟨ 0 2 0 0 0 0 5 ⟩ 5  
╵ 0 0 1 1   ╵ 0 0 1 2   ╵ 0 0 2 2                        
  0 1 0 0     0 3 0 0     0 7 0 0                        
  0 1 1 0     0 4 5 0     0 7 7 0                        
  1 1 0 0     6 7 0 0     7 7 0 0                        
          ┘           ┘           ┘                      

Python 3, 120 bytes

I,i,*g=input(),0
for c in I:g+={*sum(([*k]for k in g if{i-I.find(','),i}&k),[i+1])*(','>c)},;i+=1
print(max(map(len,g)))

Try it online!

Port of my answer to Count Islands.

Takes a comma-separated string with . representing water and * representing land as input.

-1 byte thanks to Lucenaposition

J, 42 bytes

0>./@,1#.[:+./ .*^:_~1>:[:|@-/~$j./@#:I.@,

Try it online!

• Coordinates of 1s as complex numbers
• Table of distances from each other
• Where distance is <= 1
• Transitive closure of that matrix under self multiplication
• Sum of rows
• Max

R, 108 102 bytes

-6 thanks to @pajonk

\(x)max(`[<-`(table(Reduce(\(y,i)t(apply(pmax(y,rbind(y[-1,],0))*!!y,2,rev)),1:1e3,x*seq(!x))),"0",0))

Attempt This Online!

Unpacked:

\(x)
  max( # chose biggest island
    `[<-`( # change # of 0's to 0
      table( # count island sizes
        Reduce(
          \(y,i)
            t(apply( # rotate 90 degrees
              pmax(y, rbind(y[-1,], 0)) * # if the tile directly below has bigger value, replace this tile with that tile
                (y>0), # but only if this tile isn't 0
            2, rev)),
          1:1e3,
          x * 1:length(x) # init - x, with land replaced with unique #s
        )
      ),
    "0", 0)
  )

MATL, 15 bytes

tz?4&1ZIXz&XM}z

Try it at MATL online! Or verify all test cases.

Explanation

t         % Implicit input. Duplicate
z         % Number of non-zeros
?         % If not zero
  4&1ZI   %   Label connected components with distinct positive integers, 
          %   using 4-neighbourhood to define connectedness
  Xz      %   Column vector of all the nonzero values
  &XM     %   Second output of the 'mode' function: gives the number of
          %   occurrences of the mode
}         % Else (the input is empty or only contains zeros)
  z       %   Number of non-zeros: gives 0
          % Implicit end. Implicit display

Uiua, 7 bytes

⬚0/↥⊜⧻.

2d ⊜ my beloved <3

Try it!

Get the ⧻ length of each ⊜ partition (contiguous area >0), and / reduce with ↥ maximum, defaulting to 0.

Maple, 148 bytes

proc(M)s:=[lhs]~(op(2,M));max(0,map2(op,1,LinearAlgebra:-StronglyConnectedBlocks(Matrix(nops(s)$2,(i,j)->`if`({abs~(s[i]-s[j])[]}={0,1},1,0)))))end;

Above is updated shorter solution

Input is Matrix. Finds set s of [x,y] coordinates of cells with nonzero entries (size n). Uses adjacency procedure as below to construct the nxn adjacency matrix (=1 for (i,j) adjacent, =0 otherwise). StronglyConnectedBlocks splits this into a list of submatrices that are the connected components, with no adjacencies between components. Then we find the maximum dimension of the submatrices.

Earlier approach - 177 bytes:

proc(M)                         # input is Matrix
s:=[lhs]~(op(2,M));             # set of [x,y] coords of nonzero entries
L:=table(s=~s);                 # initialize table of connectivities for each cell 
# next line iterates through, adding ones that are adjacent to any of the previously found ones
seq([for i,e in eval(L)do seq(`if`(ormap((b,c)->{abs~(b-c)[]}={0,1},[e],j),(L[i],=j),0),j=s)od],k=s);
max(0,seq(nops({i}),i=eval(L))) # count unique ones and find the max
end;

Find coordinates of cells, and initializes a table with entries for each cell. This then accumulates all other cells in the same island as that cell. Adjacency is detected by (b,c)->{abs~(b-c)[]}={0,1}. Completes by finding maximum of island sizes for each cell.

Jelly, 17 bytes

ŒṪWạ§ỊẸʋƇ@ƬṪLʋ€`Ṁ

A monadic Link that accepts a list of lists of 0 (sea) and 1 (land) and yields the largest island size.

Try it online! Or see the test-suite.

How?

ŒṪWạ§ỊẸʋƇ@ƬṪLʋ€`Ṁ - Link: list of lists of 0s and 1s, Map
ŒṪ                - truthy coordinates -> AllLand = [[Row, Column] for each 1 in Map]
               `  - use as both arguments of:
              €   -   for each {[Row, Column]}:
             ʋ    -     last four links as a dyad - f([Row, Column], AllLand):
  W               -         wrap -> FoundSoFar = [[Row, Column]]
          Ƭ       -         collect up until a fixed point under:
         @        -             with swapped arguments:
        Ƈ         -                 filter keep those of {AllLand} for which:
       ʋ          -                     last four links as a dyad - f(AllLand, FoundSoFar)
   ạ              -                         absolute difference (vectorises)
    §             -                         sums -> Manhatten distances
     Ị            -                         insignificant? (vectorises) -> IsSameOrNeighbour
      Ẹ           -                         any?
           Ṫ      -         tail -> all found orthogonally connected land locations
            L     -         length -> size of island that [Row, Column] is part of
                Ṁ - maximum

Charcoal, 42 bytes

Ｅθ⭆ι∧λψ⊞υＬＫＡＦＬθＦＬ⌈θ«Ｊκι¤0⊞υ⁻ＬＫＡΣυ»⎚Ｉ∨⌈Φυκ⁰

Try it online! Link is to verbose version of code. Explanation:

Ｅθ⭆ι∧λψ

Draw the input array to the canvas, but replace 1s with null bytes which can be flood filled later.

⊞υＬＫＡ

Get the initial number of 0s.

ＦＬθＦＬ⌈θ«

Loop through all potential island positions.

Ｊκι¤0

Try to flood the island with 0s.

⊞υ⁻ＬＫＡΣυ

See how much more land (if any) was flooded since last time.

»⎚Ｉ∨⌈Φυκ⁰

Delete the 0s and output the maximum amount of land flooded in one go, or 0 if the input was empty.

JavaScript (Node.js), 107 bytes

g=(s,X=o=0,Y)=>s.map((r,y)=>r.map((c,x)=>c<1?(x-X)**2+(Y-y)**2-1||g(s,x,y,k=1/Y?k+1:1,r[x]=o=o>k?o:k):0))|o

Try it online!

Port of my answer from Find the biggest chunk

05AB1E, 39 38 bytes

˜ƶsgäΔ4F¬ašøí}2Fø€ü3}*εε˜ŽqÆSèà]˜0KD¢à

Basic flood-fill algorithm I've used in a bunch of other challenges.

Try it online or verify all test cases.

Explanation:

Step 1: Convert all 1s to unique numbers:

˜          # Flatten the (implicit) input-matrix
 ƶ         # Multiply each value by its 1-based index
  s        # Swap to get the (implicit) input-matrix again
   g       # Pop and push its length (the amount of rows)
    ä      # Split the list into that many equal-sized rows

Try just step 1 online.

Step 2a: Flood-fill this matrix based on the input-matrix:

Δ         # Loop until the matrix no longer changes:

Step 2b: Add a border of 0s to the matrix:

 4F        #  Loop 4 times:
   ¬       #   Push the first row (without popping the matrix)
    a      #   Convert all values to 0 using an "isLetter" check
     š     #   Prepend that list of 0s to the matrix
      øí   #   Rotate the matrix once clockwise:
      ø    #    Zip/transpose; swapping rows/columns
       í   #    Reverse each inner row
  }        #  Close the loop

Step 2c: Convert the matrix into overlapping 3x3 blocks:

 2F        #  Loop 2 times:
   ø       #   Zip/transpose; swapping rows/columns
    €      #   Map over each inner list:
     ü3    #    Convert it into overlapping triplets
  }        #  Close the loop

Step 2d: Transform some 3x3 blocks back to 0s, based on the 0s in the input-matrix

 *         #  Multiple the 3x3 blocks to the values at the same positions of the
           #  (implicit) input-matrix

Step 2e: For each 3x3 block, get the maximum of its center and its horizontal/vertical neighbors:

 εε        #  Nested map over each 3x3 block:
   ˜       #    Flatten the 3x3 block to a list of 9 values
    ŽqÆ    #    Push compressed integer 13457
       S   #    Convert it to a list of digits: [1,3,4,5,7]
        è  #    Get the values at those indices
         à #    Pop and push the maximum
  ]        # Close the nested map, as well as the changes-loop

See this 05AB1E tip of mine (section How to compress large integers?) to understand why ŽqÆ is 13457.

Try just the first two steps online.

Step 3: Check which island is the largest, and output its size:

˜          # Flatten the flood-filled matrix
 0K        # Remove all 0s
   D       # Duplicate the list of positive integers
    ¢      # Pop both, and get the count of each
     à     # Pop and leave the largest count,
           # which is the size of the largest island
           # (which is output implicitly as result)