Skip to content

Instantly share code, notes, and snippets.

@planemad

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

Last active February 23, 2026 13:34
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save planemad/5b4fcd12e0321cb0175bf8e1ba8e1ddd to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save planemad/5b4fcd12e0321cb0175bf8e1ba8e1ddd to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. planemad revised this gist Feb 23, 2026. 2 changed files with 2346 additions and 0 deletions.
    2,167 changes: 2,167 additions & 0 deletions czmp-goa-4k-grid.geojson
    View file
    2,167 additions, 0 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
    179 changes: 179 additions & 0 deletions qgis-action-grid-creator.py
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,179 @@
    # 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)")
  2. planemad created this gist Feb 23, 2026.
    9 changes: 9 additions & 0 deletions czmp-goa-4k-grid-sample.geojson
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,9 @@
    {
    "type": "FeatureCollection",
    "name": "czmp-goa-4k-grid-sample",
    "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
    "features": [
    { "type": "Feature", "properties": { "row": 3, "col": 8 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 73.818979145225867, 15.773795299556086 ], [ 73.839516186943612, 15.773907239491878 ], [ 73.839622550312541, 15.754020262840418 ], [ 73.819089719867975, 15.753910860312464 ], [ 73.818979145225867, 15.773795299556086 ] ] ] } },
    { "type": "Feature", "properties": { "row": 26, "col": 17 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 74.006019016288249, 15.317330607016693 ], [ 74.026509711522749, 15.317415387929492 ], [ 74.026597438934246, 15.2975217558584 ], [ 74.006108824151383, 15.297439527598845 ], [ 74.006019016288249, 15.317330607016693 ] ] ] } }
    ]
    }