I am trying to tile a polygon with a minimal number of fixed sized squares.

Currently I am creating a fishnet over my polygon then spatially joining the squares that intersect the polygon. This is not optimal.

Also note the squares can be shifted vertically/horizontally but not rotated.

My end goal: The polygon represents a clipped image and I want a tiling of the clipped image. Where each tile is 300px by 300px. Some overlap is fine if that makes the problem easier.


enter image description here

Manually optimized

enter image description here

Are there any tools or algorithms that would help with this? I am proficient with Python and ArcMap.

The polygon was created from line features so I also have access to those segments. Maybe that would will help for generating the squares.

The entire polygon must be covered.

  • Would selecting only the squares where your polygon intersects its centroid be appropriate? This could significantly reduce your square counts. – Trevor J. Smith Oct 2 '18 at 17:27
  • I need to cover all of the polygon and I think that may leave uncovered areas on the outer sides of the polygon. I will update my question to include this. – Squanchy Oct 2 '18 at 17:38

I think the class of problem that you are looking at is called Polygon covering:

A covering of a polygon is a set of primitive units (e.g. squares) whose union equals the polygon. A polygon covering problem is a problem of finding a covering with a smallest number of units for a given polygon. This is an important class of problems in computational geometry. There are many different polygon covering problems, depending on the type of polygon being covered and on the types of units allowed in the covering.

While the "near fishnet" that you are trying to create would be useful, as an example, for GIS professionals trying to create indexes for map series which minimize the number of pages while maximizing the feature area per page, I think that you may find more potential answerers by tailoring this for a different audience and posting it at either of the maths sites where I have seen at least two similar sounding questions:

| improve this answer | |

In your manual optimalisation you shift your neighbouring tiles north-south and east-west in a manner where the resulting position does not lead to overlaps. As @PolyGeo says this is part of a broader mathematical challenge.

However, this is also a practical challenge for us dealing with GIS. Dealing with it from a more practical approach is therefore something we can do without dealing with the mathematical side of it.

In your example you are dealing with grids intersecting a polygon stretching out like a river network or similar. Using a fishnet will give you a less than ideal coverage. You solved this by shifting them both north-south and east-west. In a non-branching polygon, or where your branches meet, you may not have that option. At least not both of them at the same time. The below example is from a plate plan for making maps in Kenya. In QGIS this is called atlases and in ArcGIS Pro it is called data driven pages.

Kenya plate plan

If your tiles are tight like in the above example you will have the option of shifting horisontal or vertical lines of tiles. Your choices will have to be based on to what extent, or percentage, your tiles cover your original polygon.

Having made the fishnet procedure and then selected those tiles intersecting with your polygon you will have to continue based on the following alternatives - or a combination of them:

  1. Shift the tiles in each lines between extremes (north/south or east/west) and calculate which shifts will give you fewer tiles.
  2. Consider changing the tile sizes so that you can maintain the full coverage. If you are looking at producing maps based on the tiles, some overlap is usually not a problem. Neither would a proportional resizing of your tiles be a problem.

Programatically it is easy to create the fishnet and then select those overlapping with your polygon. I did this in FME from Safe Software and it was relatively straightforward. This can also be done in ArcGIS and QGIS.

Creating iterations which support the procedures 1 and 2 above is more demanding. But when dealing with GIS your "solution space" is probably broader than what would be the case in the ideal focused mathematical approaches mentioned by @PolyGeo.

A last option to consider is to rotate the grid or individual tiles to get an optimal coverage. For map purposes this can be done to some degree. Calculating this for the whole tileset would mean setting up an iteration which calculates the above 1 and/or 2 based on grid tilt. Individual rotation could also be calculated.

All the above calculations could be subject to brute force calculations meaning you just iterate through your solution options. But remember each new "dimension" to your approaches (shifting/tile resizing/rotations) increases your calculation effort exponentially.

| improve this answer | |

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.