4

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

enter image description here

4
  • It can't be rotated? What did you try?
    – BERA
    May 18, 2021 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 :)
    – sailestim
    May 18, 2021 at 13:20
  • Updated solution to get the square out of the rectangle.
    – Babel
    May 18, 2021 at 19:52
  • Also related topic gis.stackexchange.com/a/327005/99589
    – Taras
    May 19, 2021 at 9:23

1 Answer 1

7

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.

enter image description here

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: enter image description here

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:

enter image description here

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:

enter image description here

  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

enter image description here

  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: enter image description here

  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.

    enter image description here

  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: enter image description here

  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,
                        radians (45)
                    )
                ),
                boundary ($geometry)
            ),
            centroid($geometry)
        ),
        extend (
            @diagonal_half,
            0,
            length (@diagonal_half)
        )
    ),
    make_square (
        start_point (@diagonal),
        end_point (@diagonal)
    )
)
5
  • 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?
    – sailestim
    May 18, 2021 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
    – Babel
    May 19, 2021 at 5:26
  • See here a step by step guide how to build a more complex expression: gis.stackexchange.com/a/380390/88814
    – Babel
    May 19, 2021 at 5:42
  • 3
    Looks like a part of dissertation 😊
    – Taras
    May 19, 2021 at 8:55
  • 2
    If someone's willig to give me a title for it: why not ;-)
    – Babel
    May 19, 2021 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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