Yes I used a polygon library so what... I did eyeball the shape and guessed what the result should have been like but still more pleasing to write something that applies generally. Although I did put an area threshold (eyeballed) which subsets the corners to test, it only reduces run time to about 1/2 (from ~10s to ~5s) so that is not necessarily very critical. If I hadn't used this library, I would probably do sth along the lines of defining my own rectangle classes with unions, intersections etc so wouldn't be a more clever approach but much more time consuming.

from pathlib import Path

import numpy as np
import shapely

cwd = Path(__file__).parent

def parse_input(file_path):
  with file_path.open("r") as fp:
    data = list(map(lambda x: list(map(int,x.split(','))), fp.readlines()))

  return np.array(data)


def construct_shapes(coordinates, threshold):

  Itriu = np.triu_indices(coordinates.shape[0], k=2)
  squares = []

  for i0,i1 in zip(*Itriu):

    c0 = tuple(coordinates[i0,:])
    c1 = tuple(coordinates[i1,:])
    area = np.prod(abs(np.array(c0) - np.array(c1)) + np.array([1,1]))

    if area>threshold:
      c2 = (c0[0],c1[1])
      c3 = (c1[0],c0[1])
      squares.append(shapely.Polygon((c0,c3,c1,c2)))

  polygon = shapely.Polygon(coordinates)

  return polygon, squares


def solve_problem(file_name, redgreen=False, threshold=0):

  coordinates = parse_input(Path(cwd, file_name))

  if not redgreen:
    areas = np.prod(abs(coordinates[None,:] - coordinates[:,None]) +\
                    np.array([1,1])[None,None,:], axis=-1)
    max_area = np.max(areas)

  else:
    polygon, squares = construct_shapes(coordinates, threshold)
    max_area = -np.inf

    for inds,square in enumerate(squares):
      if square.area==0:
        continue

      if polygon.contains(square):
        c = np.array(list(zip(*square.exterior.coords.xy)))
        if (a:=np.prod(abs(c[0] - c[2]) + np.array([1,1])))>max_area:
          max_area = a

  return int(max_area)



if __name__ == "__main__":

  assert solve_problem("test_input") == 50
  assert solve_problem("input") == 4763509452
  assert solve_problem("test_input", True, 1) == 24
  assert solve_problem("input", True, 15000*80000) == 1516897893 # threshold by eyeballing the shape