Select Page

GIS Copilot Demo (V1.0): Random elevation and nearest neighbours

Task

Create a histogram to show elevation of 100 points, and also connect each point to its nearest 5 points.

Dataset

DEM of Pennsylvania.

Result

Geoprocessing workflow

 

Generated code:

  1. # Purpose:
  2. # – Generate 100 random points over the DEM extent (EPSG:32128), sample elevation from DEM, and store as a shapefile
  3. # – Compute 5 nearest neighbors for each point using qgis:distancematrix
  4. # – Build origin–neighbor connection lines (one line per pair) using native:pointstopath
  5. # – Create an elevation histogram (HTML) for the 100 points using qgis:vectorlayerhistogram
  6. # – All outputs are saved to the specified workspace, loaded into QGIS (except the HTML, whose path is printed)
  7. #
  8. # Notes:
  9. # – Distances are computed in a projected CRS (EPSG:32128), so no reprojection is required
  10. # – Implements stepwise processing with short waits to ensure files are ready before the next step
  11. # – Uses selected tools: qgis:distancematrix, native:pointstopath, qgis:vectorlayerhistogram
  12. # – DEM sampling is performed with rasterio
  13. # – Field names are kept <= 10 characters
  14. #
  15. # Inputs:
  16. # – DEM: C:/Users/17036/Downloads/PA_DEM.tif
  17. # – Output directory:
  18. #   C:\Users\17036\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\SpatialAnalysisAgent-master\Default_workspace
  20. import os
  21. import time
  22. import random
  23. import numpy as np
  24. import pandas as pd
  25. import geopandas as gpd
  26. from shapely.geometry import Point
  27. import rasterio
  28. from qgis.core import QgsVectorLayer, QgsProject
  29. import processing
  31. def create_elevation_knn_visualization():
  32.     # Paths
  33.    dem_path = r“C:/Users/17036/Downloads/PA_DEM.tif”
  34.     out_dir = r“C:\Users\17036\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\SpatialAnalysisAgent-master\Default_workspace”
  36.     # Helpers
  37.    def unique_path(directory, base_name, ext):
  38.         base = os.path.join(directory, f“{base_name}{ext}”)
  39.         if not os.path.exists(base):
  40.             return base
  41.         i = 1
  42.         while True:
  43.             cand = os.path.join(directory, f“{base_name}_{i}{ext}”)
  44.             if not os.path.exists(cand):
  45.                 return cand
  46.             i += 1
  48.     # 1) Create 100 random points within DEM bounds and sample elevation
  49.    with rasterio.open(dem_path) as src:
  50.         bounds = src.bounds  # (left, bottom, right, top)
  51.        dem_crs_epsg = 32128  # From provided overview
  52.        rng = np.random.default_rng(42)
  53.         xs = rng.uniform(bounds.left, bounds.right, 100)
  54.         ys = rng.uniform(bounds.bottom, bounds.top, 100)
  55.         coords = list(zip(xs, ys))
  56.         # Sample the DEM at point coordinates
  57.        elev_vals = list(src.sample(coords))
  58.         elev_vals = [float(v[0]) if v.size > 0 else np.nan for v in elev_vals]
  60.     points = [Point(xy) for xy in coords]
  61.     df = pd.DataFrame({
  62.         “id”: np.arange(1, 101, dtype=int),
  63.         “elev”: elev_vals
  64.     })
  65.     gdf = gpd.GeoDataFrame(df, geometry=points, crs=f“EPSG:{dem_crs_epsg}”)
  67.     # Save points shapefile with elevation
  68.    points_path = unique_path(out_dir, “points_100”, “.shp”)
  69.     gdf.to_file(points_path, driver=“ESRI Shapefile”)
  71.     # Load points into QGIS
  72.    points_vl = QgsVectorLayer(points_path, “Points_100”, “ogr”)
  73.     QgsProject.instance().addMapLayer(points_vl)
  74.     time.sleep(1.5)
  76.     # 2) Reproject for distance calculations – Not needed (already in projected CRS EPSG:32128)
  78.     # 3) k-Nearest Neighbors: compute 5 nearest neighbors using qgis:distancematrix
  79.    # Ensure we use a valid field name from layer
  80.    pts_fields = [f.name() for f in points_vl.fields()]
  81.     input_id_field = “id” if “id” in pts_fields else pts_fields[0]
  83.     dist_matrix_path = unique_path(out_dir, “knn5_matrix”, “.shp”)
  84.     dm_params = {
  85.         “INPUT”: points_vl,
  86.         “INPUT_FIELD”: input_id_field,
  87.         “TARGET”: points_vl,
  88.         “TARGET_FIELD”: input_id_field,
  89.         “MATRIX_TYPE”: 0,  # Linear (N*k x 3)
  90.        “NEAREST_POINTS”: 5,
  91.         “OUTPUT”: dist_matrix_path
  92.     }
  93.     dm_result = processing.run(“qgis:distancematrix”, dm_params)
  94.     dm_layer = QgsVectorLayer(dm_result[“OUTPUT”], “KNN5_Matrix”, “ogr”)
  95.     QgsProject.instance().addMapLayer(dm_layer)
  96.     time.sleep(1.5)
  98.     # 4) Create connection lines from origin–neighbor pairs using native:pointstopath
  99.    #    We will build a 2-point-per-pair point layer with fields: pairid (int), ord (int), dist (float)
  100.    #    GROUP_EXPRESSION = ‘pairid’, ORDER_EXPRESSION = ‘ord’
  101.    # Read back distance matrix rows
  102.    dm_fields = [f.name() for f in dm_layer.fields()]
  103.     # Find likely field names (case-insensitive fallback)
  104.    def find_field(candidates, fallback=None):
  105.         for cand in candidates:
  106.             for name in dm_fields:
  107.                 if name.lower() == cand.lower():
  108.                     return name
  109.         return fallback
  111.     in_field = find_field([“InputID”, “Input”, “INPUTID”, “INPUT”, “in_id”, “id1”], None)
  112.     tg_field = find_field([“TargetID”, “Target”, “TARGETID”, “TARGET”, “tg_id”, “id2”], None)
  113.     dist_field = find_field([“Distance”, “DISTANCE”, “dist”, “Dist”], None)
  115.     # Build id -> shapely Point mapping from gdf
  116.    id_to_geom = dict(zip(gdf[“id”].astype(int).tolist(), gdf.geometry.tolist()))
  118.     pair_rows = []
  119.     pair_counter = 1
  120.     for feat in dm_layer.getFeatures():
  121.         try:
  122.             oid = int(feat[in_field])
  123.             tid = int(feat[tg_field])
  124.             dval = float(feat[dist_field])
  125.             # Origin and target geometries
  126.            p_o = id_to_geom.get(oid, None)
  127.             p_t = id_to_geom.get(tid, None)
  128.             if (p_o is None) or (p_t is None):
  129.                 continue
  130.             # Two rows for this pair
  131.            pair_rows.append({“pairid”: pair_counter, “ord”: 0, “dist”: dval, “o_id”: oid, “t_id”: tid, “geometry”: p_o})
  132.             pair_rows.append({“pairid”: pair_counter, “ord”: 1, “dist”: dval, “o_id”: oid, “t_id”: tid, “geometry”: p_t})
  133.             pair_counter += 1
  134.         except:
  135.             # Skip problematic rows silently
  136.            continue
  138.     pairs_gdf = gpd.GeoDataFrame(pair_rows, geometry=“geometry”, crs=gdf.crs)
  139.     pairs_points_path = unique_path(out_dir, “pairs_pts”, “.shp”)
  140.     pairs_gdf.to_file(pairs_points_path, driver=“ESRI Shapefile”)
  141.     pairs_vl = QgsVectorLayer(pairs_points_path, “Pairs_Points”, “ogr”)
  142.     QgsProject.instance().addMapLayer(pairs_vl)
  143.     time.sleep(1.5)
  145.     # Run native:pointstopath to create a line per pairid
  146.    lines_path = unique_path(out_dir, “knn5_lines”, “.shp”)
  147.     p2l_params = {
  148.         “INPUT”: pairs_points_path,
  149.         “CLOSE_PATH”: False,
  150.         “ORDER_EXPRESSION”: “ord”,
  151.         “NATURAL_SORT”: False,
  152.         “GROUP_EXPRESSION”: “pairid”,
  153.         “OUTPUT”: lines_path,
  154.         “OUTPUT_TEXT_DIR”: None
  155.     }
  156.     lines_result = processing.run(“native:pointstopath”, p2l_params)
  157.     lines_vl = QgsVectorLayer(lines_result[“OUTPUT”], “KNN5_Lines”, “ogr”)
  158.     QgsProject.instance().addMapLayer(lines_vl)
  159.     time.sleep(1.5)
  161.     # 5) Histogram: Generate a histogram of the elevation attribute for the 100 points
  162.    hist_out = unique_path(out_dir, “elev_histogram”, “.html”)
  163.     # Ensure the elevation field name is correct
  164.    pts_fields = [f.name() for f in points_vl.fields()]
  165.     elev_field = “elev” if “elev” in pts_fields else pts_fields[0]
  166.     hist_params = {
  167.         “INPUT”: points_vl,
  168.         “FIELD”: elev_field,
  169.         “BINS”: 12,
  170.         “OUTPUT”: hist_out
  171.     }
  172.     hist_result = processing.run(“qgis:vectorlayerhistogram”, hist_params)
  174.     # Print only the HTML file path per requirements
  175.    print(hist_result[“OUTPUT”])
  178. # Execute
  179. create_elevation_knn_visualization()