I had a go at this last night, with this polygon:

The code I used added all the rectangles created in the process (drawn using a color based on the area):

The largest polygon was this on:

Which doesn't look that bad, although there are many things that will go wrong with the code. The idea I used was to get a first point on the boundary of the polygon, and select a second point on the boundary and have the points move around to get many combinations. These two points are treated as a side of the rectangle to be generated (one should also consider that the points could be used as diagonal of the rectangle). For this side, the perpendicular lines at start and end are created and intersected with the polygon (you should consider multipart results). The shortest line length is applied to the other line and the rectangle is created.

This is the code I came up with:

```
def main():
import arcpy
arcpy.env.overwriteOutput = True
fc = r'C:\GeoNet\MaximumAreaInscribedRectangle\data.gdb\polygon'
fc_out = r'C:\GeoNet\MaximumAreaInscribedRectangle\data.gdb\rectangles01'
polygon = arcpy.da.SearchCursor(fc, ('SHAPE@')).next()[0]
sr = polygon.spatialReference
polyline = polygon.boundary()
extent = polygon.extent
diagonal_length = getPythagoras(extent.width, extent.height)
lst_rectangles = []
for i in range(10):
perc_p1 = float(i) / 10
pntg1 = polyline.positionAlongLine(perc_p1, True)
for j in range(10):
frac = float(j) / 20
perc_p2 = perc_p1 + 0.25 + frac
if perc_p2 > 1:
perc_p2 -= 1
pntg2 = polyline.positionAlongLine(perc_p2, True)
angle = getAngle(pntg1, pntg2)
pntg1a = pntg1.pointFromAngleAndDistance(angle + 90, diagonal_length, 'PLANAR')
pntg1b = pntg1.pointFromAngleAndDistance(angle - 90, diagonal_length, 'PLANAR')
line1 = createStraightLine(pntg1a, pntg1b)
pntg2a = pntg2.pointFromAngleAndDistance(angle + 90, diagonal_length, 'PLANAR')
pntg2b = pntg2.pointFromAngleAndDistance(angle - 90, diagonal_length, 'PLANAR')
line2 = createStraightLine(pntg2a, pntg2b)
line1cut = checkInvertedLine(line1.intersect(polygon, 2), pntg1)
line2cut = checkInvertedLine(line2.intersect(polygon, 2), pntg2)
length1 = line1cut.length
length2 = line2cut.length
if length2 < length1:
line1cut = line1cut.segmentAlongLine(0, length2, False)
else:
line2cut = line2cut.segmentAlongLine(0, length1, False)
rectangle = createRectanglePolygon(line1cut, line2cut)
lst_rectangles.append(rectangle)
arcpy.CopyFeatures_management(lst_rectangles, fc_out)
def createRectanglePolygon(line1, line2):
sr = line1.spatialReference
pnt1a = line1.firstPoint
pnt1b = line1.lastPoint
pnt2a = line2.firstPoint
pnt2b = line2.lastPoint
return arcpy.Polygon(arcpy.Array([pnt1a, pnt1b,
pnt2b, pnt2a,
pnt1a]), sr)
def checkInvertedLine(line, pntg):
sr = line.spatialReference
pnt_start = line.firstPoint
pnt_end = line.lastPoint
dist_start = getPythagoras(pnt_start.X-pntg.firstPoint.X,
pnt_start.Y-pntg.firstPoint.Y)
dist_end = getPythagoras(pnt_end.X-pntg.firstPoint.X,
pnt_end.Y-pntg.firstPoint.Y)
if dist_end < dist_start:
return arcpy.Polyline(arcpy.Array([pnt_end, pnt_start]), sr)
else:
return line
def createStraightLine(pntg1, pntg2):
sr = pntg1.spatialReference
return arcpy.Polyline(arcpy.Array([pntg1.firstPoint, pntg2.firstPoint]), sr)
def getAngle(pntg1, pntg2):
'''Define angle between two points'''
return pntg1.angleAndDistanceTo(pntg2, 'PLANAR')[0]
def getPythagoras(w, h):
import math
return math.hypot(w, h)
if __name__ == '__main__':
main()
```

It is not much, but maybe a start...