Create an inscribed square inside a polygon

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 • It can't be rotated? What did you try?
– BERA
May 18 '21 at 13:06
• @BERA It can't be rotated because of the following data manipulations. Actually I don't have any ideas about the solution. I had one about the fishnet of points inside the polygon and the center of mass, but don't think it's good :) May 18 '21 at 13:20
• Updated solution to get the square out of the rectangle. May 18 '21 at 19:52
• Also related topic gis.stackexchange.com/a/327005/99589 May 19 '21 at 9:23

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:

1. Use Geometry by expression to create a bounding box, subdivided in 4 sections centered on the pole of inaccessibility with this expression:
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: 1. Menu Vector / Geometry Tools / Multipart to singleparts

2. Menu Vector / Geoprocessing Tools / Intersection, set the polygon as Input layer and the singleparts from step 2 as Overlay layer 1. Menu Vector / Geometry tools / Polygons to lines

2. Add a unique id to the layer create in step 4, using field calculator with \$id or @row_number to create a field id

3. Menu procesing / Toolbox / Explode lines

4. 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))))

Result: the diagonals are selected: 5. Toggle editing, delete selected (diagonals)

6. 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. 7. Use Geometry by expression to create the pole of inaccessibility of the original polygon with this expression: pole_of_inaccessibility (\$geometry,1)

8. 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: 1. Use Geometry by expression with the following expression: bounds(\$geometry): you get the rectangle from the very first screenshot above.

2. 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,
)
),
boundary (\$geometry)
),
centroid(\$geometry)
),
extend (
@diagonal_half,
0,
length (@diagonal_half)
)
),
make_square (
start_point (@diagonal),
end_point (@diagonal)
)
)
• Wow! That's great, thanks! I am really surprised with your elegant and complex solution in one. Could you please suggest any courses or literature for the young GIS analysts to improve the level of work with expressions and geometry? May 18 '21 at 22:02
• I learned using expressions reading the help in the expression string builder - just select a function in the middle part to see a description on the right side. Than try building an expression step by step. Also this site here helped me a lot when I was stucked. For a tutorial, see youtu.be/IXPCec8vgLA May 19 '21 at 5:26
• See here a step by step guide how to build a more complex expression: gis.stackexchange.com/a/380390/88814 May 19 '21 at 5:42
• Looks like a part of dissertation 😊 May 19 '21 at 8:55
• If someone's willig to give me a title for it: why not ;-) May 19 '21 at 9:02