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

Its classic task for https://en.m.wikipedia.org/wiki/Bisection_method — FelixIP, Aug 09 '20 at 03:42
I used Python and arcgis buffer to solve it. I can find my solution if you are interested in my answers. But language doesn't matter, algorithm is universal and very simple and very fast. — FelixIP, Aug 09 '20 at 03:53
https://gis.stackexchange.com/questions/215813/increase-each-polygon-area-to-the-same-area — FelixIP, Aug 09 '20 at 03:57
hmmm, arcpy code seems complex than i expected, i donot have arcgis license — Roman, Aug 09 '20 at 04:14

Leo · Accepted Answer · 2020-08-09T10:09:27.757

3

The .area in shapely is plane area, so we need to transform longitude, latitude (4326) to project spatial reference (e.g. 3857).
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

edited Aug 09 '20 at 10:09

answered Aug 09 '20 at 03:08

Leo

66
3

I tried according to you and edited the question. Unfortunately, I did not get expected result – Roman Aug 09 '20 at 03:22
Sorry. I think you can project the polygon first and then calculate the area. – Leo Aug 09 '20 at 03:24
how can i project? – Roman Aug 09 '20 at 03:27
Both polygons are in latlon coordinates – Roman Aug 09 '20 at 03:33
And both polygons are rectangles, so a different trick, rather than .buffer, should work, as output is enlarged polygon – Roman Aug 09 '20 at 04:03
I update the python scripts. This is for rectangles. Note: area in shapely is plane area. – Leo Aug 09 '20 at 05:30
Can you explain what the equation of buffer_size is doing? – Roman Aug 09 '20 at 06:15
Updated. As the polygon is rectangle, there is no overlaps between buffers. So we can calculate the buffer distance directly. – Leo Aug 09 '20 at 06:32
Don't use epsg:3857 for measurements it distorts area by latitude! – Ian Turton Aug 09 '20 at 09:14
Correct. Update to EPSG: 8858 (WGS 84 / Equal Earth Americas) – Leo Aug 09 '20 at 10:10
@IanTurton Which epsg is better for Asia for this purpose? – Roman Aug 09 '20 at 11:31
@Luo You did not use projected coordinates in your equation, why? – Roman Aug 09 '20 at 12:08

Buffering a polygon to match an area

1 Answers1