#!/usr/bin/env python
from osgeo import ogr
from osgeo import osr
from osgeo import gdal
import os.path

t_srs = osr.SpatialReference()
t_srs.ImportFromEPSG(26915)

s_srs = osr.SpatialReference()
s_srs.ImportFromEPSG(4326)

def raster_extents(raster):
    ds=ogr.Open(raster)

    layer=ds.GetLayerByIndex(0)
    left, right, top, bottom = layer.GetExtent()

    feature = layer.GetNextFeature()
    while feature:
        geom = feature.GetGeometryRef().Clone()
        geom.TransformTo(t_srs)
        yield geom

        feature = layer.GetNextFeature()

extents = {}
import glob
for outline in glob.glob('cutlines/*.json'):
    for ext in raster_extents(outline):
        fn, _ = os.path.splitext(os.path.basename(outline))
        vrt = 'data/%s.vrt' % (fn,)
        extents[vrt] = ext

for x1 in range(460000, 520000, 10000):
    for y1 in range(4920000, 4970000, 10000):
        left = x1
        top = y1
        right = x1 + 10000
        bottom = y1 + 10000

        edge = ogr.Geometry(ogr.wkbLinearRing)
        edge.AssignSpatialReference(t_srs)
        edge.AddPoint(left, top)
        edge.AddPoint(right, top)
        edge.AddPoint(right, bottom)
        edge.AddPoint(left, bottom)
        edge.CloseRings()

        cutline = ogr.Geometry(ogr.wkbPolygon)
        cutline.AddGeometry(edge)

        entries = []
        for k, v in extents.iteritems():
            #print("Edge", str(cutline))
            #print("V", str(v))
            if not v.Disjoint(cutline):
                entries.append(k)

        print("%d %d (%d)" % (x1, y1, len(entries)))
        if len(entries):
            with open('tile-entries/%dx%d.txt' % (x1, y1), 'w') as out:
                for i in entries:
                    out.write(i + '\n')

@@ -98,16 +98,150 @@
 s.path%0A%0A
+STEP = 10000%0AMINX = 180000 #189774.764105%0AMINY = 4810000 #4816337.325688%0AMAXX = 770000 #761944.028930%0AMAXY = 5480000 #5472405.931701%0A%0A
 t_srs = 
@@ -938,29 +938,24 @@
 nge(
-460000, 520000, 10000
+MINX, MAXX, STEP
 ):%0A 
@@ -977,31 +977,24 @@
 nge(
-4920000, 4970000, 10000
+MINY, MAXY, STEP
 ):%0A 
@@ -1048,21 +1048,20 @@
  = x1 + 
-10000
+STEP
 %0A       
@@ -1075,21 +1075,20 @@
  = y1 + 
-10000
+STEP
 %0A%0A      
@@ -1630,16 +1630,45 @@
 end(k)%0A%0A
+        if len(entries):%0A    
         
@@ -1716,33 +1716,8 @@
 )))%0A
-        if len(entries):%0A
     

Factor out min/max bounds