2

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

6

1 Answer 1

3
  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



12
  • I tried according to you and edited the question. Unfortunately, I did not get expected result
    – Roman
    Aug 9, 2020 at 3:22
  • Sorry. I think you can project the polygon first and then calculate the area.
    – Leo
    Aug 9, 2020 at 3:24
  • how can i project?
    – Roman
    Aug 9, 2020 at 3:27
  • Both polygons are in latlon coordinates
    – Roman
    Aug 9, 2020 at 3:33
  • And both polygons are rectangles, so a different trick, rather than .buffer, should work, as output is enlarged polygon
    – Roman
    Aug 9, 2020 at 4:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.