Creating polygon by specifying area and some (irregular) boundaries

Question

To assist with a multi-year planting project, I need to find a way to identify and map the portion of the whole planting area that will be planted in year 1.

With reference to the image below, the things I know are:

The area we will plant in year 1 is 1.4 Ha.
The total planting area (orange polygon) for the multi year project is approximately 10 times larger.
We will start planting at A, in the north, and work our way south filling the orange polygon as we go.
We will plant out to the edges of the orange polygon identified by the solid blue lines.
What I want to know in advance is:

Where is the dashed blue line going to end up?

How much of the irregularly shaped orange polygon amounts to 1.4Ha of planting?

I want to automate this, not adjust by trial and error.

Ideally the result is a new polygon with an area of 1.4Ha that fills the northern end of my orange polygon.

Welcome to GIS SE. As a new user, please take the [Tour]. A Question like this, which is probably too open-ended for a reasonable answer, could be improved by a graphic or three, explaining what you seek to accomplish. There are obviously an infinite number of potential solutions here, and "elegant" in that context is probably too opinion-based to vote on a "best" answer. — Vince, Sep 08 '22 at 02:27
Can you add some more examples of shapes? How many polygons are there in total? Do you have a point layer with starting points for each polygon? — BERA, Sep 08 '22 at 05:41
Does this help? https://gis.stackexchange.com/a/354691/88814 — Babel, Sep 08 '22 at 06:30

score 4 · Answer 1 · edited Sep 12 '22 at 03:37

I have created a script that work well for a few examples I made.

The calculations are all cartesian, so if the Coordinate Reference System of your polygon layer has poor accuracy or is geographic, this script will be either inaccurate or won't work

Example: result of the script:

import numpy as np
Ojective to reach in meters squared
OBJECTIVE_METER_SQUARED = 14000.
Tolerance in meter squared of the objective to accept calculated result
TOLERANCE = 1.
Angle of the conic angled lines to the polygon centerline
centerline : northest point -> centroid point -> southest point
Check that the input polygon doesn't cross the two outer lines
if the input polygon crosses line, increase ANGLE_CONIC_CURVE
ANGLE_CONIC_CURVE = 30
Name of the polygon layer name in legend
POLY_NAME = "poly"
Get input polygon
poly_layer = QgsProject.instance().mapLayersByName(POLY_NAME)[0]
geom_poly = list(poly_layer.getFeatures())[0].geometry()
polygon = geom_poly.asPolygon()
Create centerline and conic curved outer lines
Find northest point and southest point
northest_point = None
southest_point = None
ymax = None
ymin = None
for list_poly_pt in polygon:
    for pt in list_poly_pt:
        if ymax is None or pt.y() > ymax:
            ymax = pt.y()
            northest_point = pt
        if ymin is None or pt.y() < ymin:
            ymin = pt.y()
            southest_point = pt
centroid_point = geom_poly.centroid().asPoint()
Create the curved center line
curve_part = QgsCircularString()
curve_part.setPoints([
    QgsPoint(northest_point),
    QgsPoint(centroid_point),
    QgsPoint(southest_point)
])
curve_feat = QgsFeature()
geom_curve = QgsGeometry(curve_part)
curve_feat.setGeometry(geom_curve)
Create the curved outer line 1
curve_feat_minus_angle = QgsFeature()
geom_minus_angle = QgsGeometry(geom_curve)
geom_minus_angle.rotate(ANGLE_CONIC_CURVE, northest_point)
curve_feat_minus_angle.setGeometry(geom_minus_angle)
Create the curved outer line 2
curve_feat_plus_angle = QgsFeature()
geom_plus_angle = QgsGeometry(geom_curve)
geom_plus_angle.rotate(-ANGLE_CONIC_CURVE, northest_point)
curve_feat_plus_angle.setGeometry(geom_plus_angle)
Export to legend the curved lines
curve_layer = QgsVectorLayer(
    "LineString?crs=epsg:4326&index=yes", 
    "temporary_curve",
    "memory"
)
curve_layer.setCrs(poly_layer.crs())
curve_layer.dataProvider().addFeatures([
    curve_feat,
    curve_feat_minus_angle,
    curve_feat_plus_angle,
])
QgsProject.instance().addMapLayer(curve_layer)
Creates a curved cone shaped polygon that follow x pourcent of the
distance of the three curved lines passed in the parameters
def get_geom_interpolated_poly(curve0, curve1, curve2, pourcent):
    dist0 = min(curve0.length()/100. * pourcent, curve0.length())
    dist1 = min(curve1.length()/100. * pourcent, curve1.length())
    dist2 = min(curve2.length()/100. * pourcent, curve2.length())
    pt_start = curve0.constGet().startPoint()
    pt0 = curve0.interpolate(dist0).asPoint()
    pt1 = curve1.interpolate(dist1).asPoint()
    pt_half1 = curve1.interpolate(dist1/2).asPoint()
    pt2 = curve2.interpolate(dist2).asPoint()
    pt_half2 = curve2.interpolate(dist2/2).asPoint()
line = QgsLineString([])
curve_part1 = QgsCircularString()
curve_part1.setPoints([
    pt_start,
    QgsPoint(pt_half1),
    QgsPoint(pt1)
])
curve_part0 = QgsCircularString()
curve_part0.setPoints([
    QgsPoint(pt1),
    QgsPoint(pt0),
    QgsPoint(pt2)
])
curve_part2 = QgsCircularString()
curve_part2.setPoints([
    QgsPoint(pt2),
    QgsPoint(pt_half2),
    pt_start,
])

curve = line.toCurveType()
curve.addCurve(curve_part1)
curve.addCurve(curve_part0)
curve.addCurve(curve_part2)
polygon = QgsPolygon()
polygon.setExteriorRing(curve)
return QgsGeometry(polygon)


def calculate_area(poly_ori, poly_cur):
    inter = poly_ori.intersection(poly_cur)
    return inter, inter.area()
pourcent_change = 10
current_area = None
current_geom = None
dif_area = None
under = (OBJECTIVE_METER_SQUARED - TOLERANCE)
over = (OBJECTIVE_METER_SQUARED + TOLERANCE)
pourcent_under = 0
pourcent_over = 100
While the current area calculated is not close enough to the objective
the pourcent_change gets smaller and smaller unter the area calulated
is under the tolerance from the objective area
while current_area is None or not (under <= current_area <= over):
    for pourcent in np.arange(pourcent_under, pourcent_over+1, pourcent_change):
        cur_geom = get_geom_interpolated_poly(
            geom_curve,
            geom_minus_angle,
            geom_plus_angle,
            pourcent
        )
        current_geom , current_area = calculate_area(geom_poly, cur_geom)
        if current_area <= OBJECTIVE_METER_SQUARED:
            pourcent_under = pourcent
        else:
            pourcent_over = pourcent
            break
    pourcent_change /= 10
Create polygon feature
poly_feat = QgsFeature()
poly_feat.setGeometry(current_geom)
Create a memory layer
layer_out = QgsVectorLayer(
    "Polygon?crs=epsg:4326&index=yes",
    "Output Polygon",
    "memory"
)
layer_out.setCrs(poly_layer.crs())
layer_out.dataProvider().addFeatures([poly_feat])
QgsProject.instance().addMapLayer(layer_out)

Final area will always be between OBJECTIVE_METER_SQUARED and OBJECTIVE_METER_SQUARE + TOLERANCE, never under (for you between 14000 and 14001 if tolerance is 1)

If you want a flat lined polygon output at the border of the polygon you can change how the polygons are created inside the get_geom_interpolated_poly function. I created the algorithm as I did because it seemed the "most accurate" to create a representation of the propagation (progress) from a single source evolving in a polygon

Incredible solution. This should become a processing tool. Maybe with the possibility to not only have one single feature as output but to divide the whole polygon into segments of same size (until the last which will take the rest). — Bernd Loigge, Sep 08 '22 at 09:54
To allow more complex shapes or a wider start of the polygon you could offset the start (northest_point) of the curved lines to allow the lines have more distance to "spread" and the polygon won't cross the lines — Kalak, Sep 08 '22 at 09:54
@BerndLoigge Could be interesting, I don't think it would even be that hard, a QgsMapTool tool to get the origin of the direction of the propagation, some arguments to decide the shape (flat slices, curved slices) of the propagation , rotation, .. Maybe another argument for deciding if the splits should be made every x % of the polygon area or x m² of the polygon area — Kalak, Sep 08 '22 at 10:42
@Alasdair your welcome. Remember to accept the answer please when a satisfactory answer has been given — Kalak, Sep 09 '22 at 05:06

Creating polygon by specifying area and some (irregular) boundaries

1 Answers1

Ojective to reach in meters squared

Tolerance in meter squared of the objective to accept calculated result

Angle of the conic angled lines to the polygon centerline

centerline : northest point -> centroid point -> southest point

Check that the input polygon doesn't cross the two outer lines

if the input polygon crosses line, increase ANGLE_CONIC_CURVE

Name of the polygon layer name in legend

Get input polygon

Create centerline and conic curved outer lines

Find northest point and southest point

Create the curved center line

Create the curved outer line 1

Create the curved outer line 2

Export to legend the curved lines

Creates a curved cone shaped polygon that follow x pourcent of the

distance of the three curved lines passed in the parameters

While the current area calculated is not close enough to the objective

the pourcent_change gets smaller and smaller unter the area calulated

is under the tolerance from the objective area

Create polygon feature

Create a memory layer