Geographic Information Systems Asked by sailestim on May 21, 2021
I have polygons of the administrative boundaries and want to create a square polygon inside of them, which will have the maximum area and will be fully covered by these boundaries. In other words: how to create an inscribed square inside a polygon.
The example is below
This solution is based on the 4 white dots in the following screenshots. They are the closest points to the or pole of inaccessability for each section of the outer boundary of the polygon. Section here means: divide the polygon in four parts (north/east/south/west) in a diagonal-like manner. From this rectangle, get the shorter side and reduce the longer side to this length.
From the four white points, create a bounding box to get the blue rectangle. From this, it's ease to get the red square, based on its diagonal: starting at the or pole of inaccessability, going in an angle of 45 degrees until it intersects the boundary of the rectangle.
Find a step-by-step description below.
Remarks:
Better than using centroid ()
(as propagated in the first version of this solution) is pole_of_inaccessibility( )
(see here for the concept of pole of inaccessability). As you can see on the following screenshot, some country shapes (like the contiguous territories of Italy or especially Croatia) are irregularily shaped in a way that the centroid does not fall in the place with most space around it inside the country borders: but we need as much space as possible to get the biggest square possible inside the territory. Thus the pole of inaccessability is better. As the first version of this solution was besed on centroid, the screenshots inserted in the steps are still based on this - with the shape used for this, centroid or pole of inaccessability are almost the same and thus does not affect much how the solution looks - for Italy or Croatia, however, it does!
Blue dots: pole of inaccessability, white dots: centroids:
Centroid of Norway even falls outside the territory of Norway: it is inside of Sweden, whereas the pole of inaccessability by definition is alwasy inside the territory:
Step by step description:
For each step, use the output of the previous step, if not stated otherwise:
collect_geometries (
array_foreach (
array (1,2,3,4),
make_triangle (
pole_of_inaccessibility ($geometry,1),
point_n( bounds( $geometry),@element),
point_n( bounds( $geometry),@element+1)
)
)
)
Orange (covered): original polygon; blue: polygons created with the expression:
Menu Vector / Geometry Tools / Multipart to singleparts
Menu Vector / Geoprocessing Tools / Intersection
, set the polygon as Input layer
and the singleparts from step 2 as Overlay layer
Menu Vector / Geometry tools / Polygons to lines
Add a unique id to the layer create in step 4, using field calculator with $id
or @row_number
to create a field id
Menu procesing / Toolbox / Explode lines
Select by expression
with this expression - replace 'poly'
with the name of your original polygon layer: disjoint ($geometry, boundary ( geometry ( get_feature_by_id ('poly',1))))
Toggle editing, delete selected (diagonals)
Menu Processing / Toolbox / Aggregate
, set the id
in the field Group by expression (NULL to group all features)
and under Source Expression
select first value
as Aggregate Function for the id
.
Use Geometry by expression to create the pole of inaccessibility of the original polygon with this expression: pole_of_inaccessibility ($geometry,1)
Use Geometry by expression on the pole of inaccessibility layer with the following expression to create the 4 closest point to the pole of inaccessibility : each for every 4 sections:
collect_geometries (
array_foreach (
array (1,2,3,4),
closest_point (
geometry (
get_feature_by_id (
'Aggregiert',
@element
)
),
$geometry
)
)
)
Centroid (or pole of inaccessability, blue) and the 4 closest point, one per each section:
Use Geometry by expression with the following expression: bounds($geometry)
: you get the rectangle from the very first screenshot above.
To convert the rectangle to a square, use Geometry by expression with the following expression - adapt the value in line 11 from 1000
in my example to your needs - it can't be too long, so if unsure, just add several zeros:
with_variable (
'diagonal',
with_variable (
'diagonal_half',
make_line (
intersection (
make_line (
centroid ($geometry),
project (
centroid ($geometry),
1000,
radians (45)
)
),
boundary ($geometry)
),
centroid($geometry)
),
extend (
@diagonal_half,
0,
length (@diagonal_half)
)
),
make_square (
start_point (@diagonal),
end_point (@diagonal)
)
)
Correct answer by Babel on May 21, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP