0

I have two overlapping pointclouds, for this example I take the first 10, A[0:10] and B[0:10]. They are stored as tuples, [(x, y, z),...].

I want to compare the two based on their nearest x,y,z neighbours, for example take A[0] and find the nearest corresponding 3D point and add the distance value to the tuple, then iterate through dataset A comparing it to Dataset B. Dataset A then has the distance value appended to the tuples [(x, y, z, dist),...].

I'm using python either externally or internally within QGIS.

Dataset A

[(580992.5136, 4275261.8512, 191.2496),
 (580992.5107, 4275261.855, 191.2295),
 (580992.5157, 4275261.8573, 191.2698),
 (580992.5208, 4275261.8582, 191.2428),
 (580992.5186, 4275261.8587, 191.2567),
 (580992.518, 4275261.8595, 191.2628),
 (580992.5115, 4275261.8597, 191.1952),
 (580992.5179, 4275261.8619, 191.2227),
 (580992.518, 4275261.8621, 191.2766),
 (580992.528, 4275261.8648, 191.2369)]

Dataset B

[(580992.4163, 4275262.2737, 191.062),
 (580992.4165, 4275262.2731, 191.0551),
 (580992.4172, 4275262.2762, 191.0681),
 (580992.4175, 4275262.2544, 191.0734),
 (580992.4196, 4275262.2696, 191.0743),
 (580992.4219, 4275262.2591, 191.0679),
 (580992.4227, 4275262.2711, 191.0618),
 (580992.4232, 4275262.2723, 191.0688),
 (580992.4249, 4275262.2711, 191.0552),
 (580992.4296, 4275262.2677, 191.0617)]

--- Edit

The points are originally stored as pgpointclouds (patches of 600 points - via PDAL), I've fetched them to QGIS and/or python using (I expect there is a better way) the following:

SELECT st_makepoint(st_x(PC_EXPLODE(pa)::geometry),st_y(PC_EXPLODE(pa)::geometry),st_z(PC_EXPLODE(pa)::geometry) ) AS geom
from pc_201807180937;
  • You say pyqgis, do you have any code other than the input data? What do you want to do with the result? Finding matches is difficult enough but what do you want to do with the matches when you have them? Export pairs to files, build a dict... – Michael Stimson Oct 1 '18 at 11:48
  • Ultimately save them in a new table in postgis which I will then filter based on the distance so I can split Dataset A into points closer to Dataset B and points far from Dataset B. Having it as a script in qgis would be very useful, but for now I'm just trying to find a solution that works. The pointclouds are fairly small ca. 5000 points each. – Spatial Digger Oct 1 '18 at 11:58
  • 1
    How about doing it in Postgis directly, using ST_ClosestPoint3D or something similar. You have the data stored in pg_pointcloud, I assume? – John Powell Oct 1 '18 at 12:04
  • John I've edited the above to answer your question. ST_ClosestPoint3D looks promising, on reading the docs it compares one line to one point, so I'm not sure on the correct syntax. I wouldn't want to manually select each starting point location, is there a looping function - this is why I was looking into Python in QGIS. – Spatial Digger Oct 1 '18 at 12:28
  • Would it be acceptable to impose some threshold distance such that points that are not within the threshold of some other point are not assigned a nearest neighbor? – Kirk Kuykendall Oct 1 '18 at 12:50
0

Updated (made some mistakes):

Out of my head, a PostGIS based solution:

WITH
    pc_a_geom AS (
        SELECT id,
               PC_PCid(pts) AS pcid,
               pts::geometry
        FROM (
            SELECT id,
                   PC_Explode(<patch_col>) AS pts
            FROM <patch_table_a>
        ) AS q
    ),

    pc_b_geom AS (
        SELECT id,
               PC_PCid(pts) AS pcid,
               pts::geometry
        FROM (
            SELECT id,
                   PC_Explode(<patch_col>) AS pts
            FROM <patch_table_b>
        ) AS q
    )

SELECT a.id,
       PC_Patch(PC_MakePoint(a.pcid, ARRAY[ST_X(a.pts), ST_Y(a.pts), ST_Z(a.pts), ST_3DDistance(a.pts, b.pts)]))
FROM pc_a_geom AS a
JOIN LATERAL (
    SELECT pts
    FROM pc_b_geom AS bb
    ORDER BY bb.pts <<->> a.pts
    LIMIT 1
) AS b
ON true
GROUP BY a.id;

with the 3D KNN operator <<->> lifting the heavy weight of the neighbor search. As there is no index, things will be slow...

Note that ST_3DDistance operates on geometry type only, thus you´d want to have a proper projection in place to get meaningful results (using EPSG:4326 will result in degrees and measurements be off between Latitudes).

This should return equal pcpatch objects as from <patch_table_a>, with all pcpoints having the distance to it´s nearest neighbor of <patch_table_b> assigned as the M value, or the 4th dimension in pgPointCloud terms.

  • ...and this takes 45 secs for a test with 5000 points in 5 patches each on a mid-tech setup. btw., make sure you have the distance dimension properly set up in the PC schema document. – ThingumaBob Oct 1 '18 at 13:50
  • Its been running for a good 20 minutes, exact numbers are patch_table_a 6132 points, patch_table_b 5115 points, stored in patches of 600 points. When you say make sure the PC schema document is setup correctly, what do you mean? I have the post_gis, pointcloud, pointcloud_postgis extensions enabled. I'm using projection espg 32635, how do I make sure I have set the distance dimension? I'm not sure what you mean by this term. – Spatial Digger Oct 1 '18 at 14:43
  • I canceled the process and added limit 100 to both of the pc_a_geom and pc_b_geom to see if it returned what I was expecting. This threw a ERROR: array dimensions do not match schema dimensions of pcid = 1 SQL state: XX000 Is this the result of adding a limit? My laptop is a new (2 weeks) MSI, i7 with 32GB Ram so it should handle this ok. – Spatial Digger Oct 1 '18 at 14:52
  • I'm now dissecting the code, the pc_a_geom and pc_b_geom parts seems fine. – Spatial Digger Oct 1 '18 at 15:06
  • the error is due to the pgpointcloud schema document that you have inserted in the pointcloud_formats table with pcid=1; all 'coordinates' of a point (called 'dimensions') are defined there and you will need to specify the data type, size and scale, especially for the new distance 'dimension' (if you have used the document from the docs, it's the intensity dimension that you have to alter) – ThingumaBob Oct 1 '18 at 15:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.