Revisions to Buffering a polygon to match an area

edited body

Source Link

edited Aug 9, 2020 at 10:09

66
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:3857'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) # 162676072602129724395664.3725697488


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

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:3857'))

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) # 162676072602.37256


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

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

added 394 characters in body

Source Link

edited Aug 9, 2020 at 6:29

Leo

66
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:3857'))

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) # 162676072602.37256


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

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:3857'))

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) # 162676072602.37256


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

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:3857'))

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) # 162676072602.37256


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

added 1046 characters in body

Source Link

edited Aug 9, 2020 at 5:26

Leo

66
3

Try: x = difference(area_poly_A - area_poly_B) / perimeter_poly_B

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:3857'))

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) # 162676072602.37256


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

Try: x = difference(area_poly_A - area_poly_B) / perimeter_poly_B

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:3857'))

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) # 162676072602.37256


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

Source Link

created Aug 9, 2020 at 3:08

Leo

66
3

Loading

Stack Exchange Network

Return to Answer