Buffering a polygon to match an area

Question

I have a large polygon (A) with a specific area. I also have another polygon (B) with area less than polygon A.

I want to buffer the polygon B so that its area becomes equal to the area of A.

from shapely.geometry import Polygon
ulx,uly= -105.645292, 28.094511
lrx, lry = -101.985052, 24.885167
poly_A = Polygon([(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)])
print (poly_A.area)


ulx,uly= -110.398408,  58.136267
lrx, lry = -108.689382,  57.137692
poly_B = Polygon([(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)])
print (poly_B.area)


buffer_size = (poly_A.area - poly_B.area) / poly_B.length
print (buffer_size)

poly_C = poly_B.buffer(buffer_size)
print (poly_C.area)

After knowing buffer_size, I can get bufferred polygon C with area equals to polygon A.

Note. All polygons dealt are rectangles. Both polygons are in latlon coordinates

Answer 1

1. The .area in shapely is plane area, so we need to transform longitude, latitude (4326) to project spatial reference (e.g. 3857).
2. For rectangle, if the original size is a * b, area is A, buffer distance is X, the buffer area B then:

B = A + 2 * (a + b) * X + pi * X * X   ----- (1)

2 * (a + b) = c = perimeter = poly_B.length

(1) => X = sqrt(c * c - 4 * (B - A) * PI) / (2 * PI)

import pyproj
import math
from functools import partial
from shapely.geometry import Polygon
from shapely.geometry import shape
from shapely.ops import transform

proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), \
               pyproj.Proj(init='epsg:8858'))

ulx, uly= -105.645292, 28.094511
lrx, lry = -101.985052, 24.885167
poly_A = Polygon([(ulx, uly),(ulx, lry),(lrx, lry),(lrx, uly)])
print (poly_A.area)  # 11.746969282560011


poly_A_proj = transform(proj, poly_A)
poly_A_proj_area = poly_A_proj.area
print(poly_A_proj_area) # 129724395664.97488


ulx, uly= -110.398408, 58.136267
lrx, lry = -108.689382, 57.137692
poly_B = Polygon([(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)])
print (poly_B.area)  # 39511382337.41971


poly_B_proj = transform(proj, poly_B)
poly_B_proj_area = poly_B_proj.area


buffer_size = (math.sqrt(poly_B.length * poly_B.length - 4*(poly_B.area - poly_A.area)*math.pi)-poly_B.length) / (2 * math.pi)
print (buffer_size)  # 1.122771717045914


poly_C = poly_B.buffer(buffer_size)
print (poly_C.area)  # 11.740610528262806