czmp-goa-4k-grid-sample.geojson

czmp-goa-4k-grid.geojson

qgis-action-grid-creator.py

# QGIS layer action Python script to construct a grid using sample frames polygons
# This can be run on a polygon layer that contains one or more identical sized grid cells for the target grid
# The output will be a complete grid with required rows x cols that matches the sample cells

from qgis.core import *
from qgis.utils import iface
import numpy as np
from PyQt5.QtCore import QVariant, Qt
from PyQt5.QtWidgets import QDialog, QFormLayout, QSpinBox, QDialogButtonBox, QLabel

# ── 1. Get all features ───────────────────────────────────────────────────────
layer = iface.activeLayer()
if layer is None:
    raise Exception("No active layer")

features = list(layer.getFeatures())
print(f"Layer: '{layer.name()}'  |  Features: {len(features)}")
if len(features) < 1:
    raise Exception("Need at least 1 feature")

# ── 2. Extract ordered corners [TL, TR, BR, BL] ──────────────────────────────
def poly_corners(feat):
    geom = feat.geometry()
    ring = geom.asMultiPolygon()[0][0] if geom.isMultipart() else geom.asPolygon()[0]
    if ring[0] == ring[-1]:
        ring = ring[:-1]
    pts = np.array([[p.x(), p.y()] for p in ring])
    cx, cy = pts[:, 0].mean(), pts[:, 1].mean()
    top = pts[pts[:, 1] >= cy]
    bot = pts[pts[:, 1] < cy]
    if len(top) < 2 or len(bot) < 2:
        xs, ys = pts[:, 0], pts[:, 1]
        return np.array([
            [xs.min(), ys.max()], [xs.max(), ys.max()],
            [xs.max(), ys.min()], [xs.min(), ys.min()]
        ])
    return np.array([
        top[top[:, 0].argmin()],  # TL
        top[top[:, 0].argmax()],  # TR
        bot[bot[:, 0].argmax()],  # BR
        bot[bot[:, 0].argmin()]   # BL
    ])

# ── 3. Step vectors from each polygon's own edges ────────────────────────────
# step_right = TR - TL  (one cell width vector)
# step_down  = BL - TL  (one cell height vector)
# Inter-frame distance is NEVER used.
all_corners = [poly_corners(f) for f in features]

step_right = np.mean([c[1] - c[0] for c in all_corners], axis=0)  # TR - TL
step_down  = np.mean([c[3] - c[0] for c in all_corners], axis=0)  # BL - TL

cell_w = np.linalg.norm(step_right)
cell_h = np.linalg.norm(step_down)

print(f"\n=== Step vectors (averaged over {len(features)} frames) ===")
print(f"  step_right : ({step_right[0]:.8f}, {step_right[1]:.8f})  |{cell_w:.6f}|")
print(f"  step_down  : ({step_down[0]:.8f}, {step_down[1]:.8f})  |{cell_h:.6f}|")

# ── 4. Back-project each frame to estimate origin (row=0, col=0 virtual TL) ──
# origin = TL - col*step_right - row*step_down
# Rows and cols are 1-based, so origin is one step up-left of cell [1,1]
origin_estimates = []
print("\n=== Origin estimates per frame ===")
for feat, corners in zip(features, all_corners):
    r = int(feat['row'])
    c = int(feat['col'])
    tl = corners[0]
    est = tl - c * step_right - r * step_down
    origin_estimates.append(est)
    print(f"  frame[r={r:3d}, c={c:3d}]  TL=({tl[0]:.6f}, {tl[1]:.6f})  "
          f"origin est=({est[0]:.6f}, {est[1]:.6f})")

origin = np.mean(origin_estimates, axis=0)
print(f"\n  Final origin: ({origin[0]:.8f}, {origin[1]:.8f})")

# ── 5. Residuals ──────────────────────────────────────────────────────────────
print("\n=== Fit residuals ===")
for feat, corners in zip(features, all_corners):
    r, c = int(feat['row']), int(feat['col'])
    predicted_tl = origin + c * step_right + r * step_down
    actual_tl    = corners[0]
    err = np.linalg.norm(predicted_tl - actual_tl)
    print(f"  frame[r={r}, c={c}]  "
          f"predicted=({predicted_tl[0]:.6f}, {predicted_tl[1]:.6f})  "
          f"actual=({actual_tl[0]:.6f}, {actual_tl[1]:.6f})  "
          f"error={err:.8f}")

all_rows = [int(f['row']) for f in features]
all_cols = [int(f['col']) for f in features]

# ── 6. Dialog ─────────────────────────────────────────────────────────────────
class GridDialog(QDialog):
    def __init__(self):
        super().__init__()
        self.setWindowTitle("Grid Builder")
        self.setWindowFlags(self.windowFlags() | Qt.WindowStaysOnTopHint)
        self.setMinimumWidth(400)
        layout = QFormLayout()

        self.start_row = QSpinBox(); self.start_row.setRange(1, 999); self.start_row.setValue(1)
        self.start_col = QSpinBox(); self.start_col.setRange(1, 999); self.start_col.setValue(1)
        self.n_rows    = QSpinBox(); self.n_rows.setRange(1, 999);    self.n_rows.setValue(max(all_rows))
        self.n_cols    = QSpinBox(); self.n_cols.setRange(1, 999);    self.n_cols.setValue(max(all_cols))

        layout.addRow("Start row (1-based):", self.start_row)
        layout.addRow("Start col (1-based):", self.start_col)
        layout.addRow("Total rows to generate:", self.n_rows)
        layout.addRow("Total cols to generate:", self.n_cols)

        layout.addRow(QLabel(
            f"\n── Input frames ──\n"
            f"  Count : {len(features)}\n"
            f"  Rows in data : {min(all_rows)} → {max(all_rows)}\n"
            f"  Cols in data : {min(all_cols)} → {max(all_cols)}\n"
            f"\n── Cell size (from polygon edges) ──\n"
            f"  Width  : {cell_w:.6f} map units\n"
            f"  Height : {cell_h:.6f} map units\n"
            f"  step_right : ({step_right[0]:.4f}, {step_right[1]:.4f})\n"
            f"  step_down  : ({step_down[0]:.4f}, {step_down[1]:.4f})\n"
            f"\n  Set total rows/cols to the FULL grid extent.\n"
            f"  e.g. 25 rows × 40 cols for 1000 sheets."
        ))

        btn = QDialogButtonBox(QDialogButtonBox.Ok | QDialogButtonBox.Cancel)
        btn.accepted.connect(self.accept)
        btn.rejected.connect(self.reject)
        layout.addRow(btn)
        self.setLayout(layout)

dlg = GridDialog()
dlg.show()
dlg.raise_()
dlg.activateWindow()

if not dlg.exec_():
    print("Cancelled.")
else:
    start_row  = dlg.start_row.value()
    start_col  = dlg.start_col.value()
    total_rows = dlg.n_rows.value()
    total_cols = dlg.n_cols.value()

    print(f"\nGenerating {total_rows} × {total_cols} = {total_rows * total_cols} cells "
          f"(rows {start_row}–{start_row + total_rows - 1}, "
          f"cols {start_col}–{start_col + total_cols - 1})")

    # ── 7. Build output layer ─────────────────────────────────────────────────
    out_layer = QgsVectorLayer(
        f"Polygon?crs={layer.crs().authid()}", "Reconstructed_Grid", "memory")
    pr = out_layer.dataProvider()
    pr.addAttributes([
        QgsField("row",    QVariant.Int),
        QgsField("col",    QVariant.Int),
        QgsField("map_id", QVariant.String),
    ])
    out_layer.updateFields()

    feats_out = []
    for r in range(start_row, start_row + total_rows):
        for c in range(start_col, start_col + total_cols):
            tl = origin + c * step_right + r * step_down
            tr = tl + step_right
            br = tr + step_down
            bl = tl + step_down
            ring = [QgsPointXY(float(p[0]), float(p[1])) for p in [tl, tr, br, bl, tl]]
            fo = QgsFeature()
            fo.setGeometry(QgsGeometry.fromPolygonXY([ring]))
            fo.setAttributes([r, c, f"R{r:03d}_C{c:03d}"])
            feats_out.append(fo)

    pr.addFeatures(feats_out)
    out_layer.updateExtents()
    QgsProject.instance().addMapLayer(out_layer)

    print(f"Done — {len(feats_out)} features in 'Reconstructed_Grid'")
    iface.messageBar().pushSuccess(
        "Grid Builder",
        f"Seamless {total_rows}×{total_cols} grid added ({len(feats_out)} cells)")