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.
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:
- Shift the tiles in each lines between extremes (north/south or east/west) and calculate which shifts will give you fewer tiles.
- 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.
