Spatial interpolation in PostGIS without outputting Raster

Question

I have some data that is similar to rainfall, ie. a latitude, longitude, and a Z number. It is in a point layer that is a nice clean point grid. Many Z numbers are missing in this grid.
I want to interpolate to find reasonable Z numbers. TIN, IDW,or other methods. I do not want to go into a Raster first then back out the intersections from the raster onto the point grid.
I can perform an interpolation in QGIS, but it generates a raster which I'm wanting to avoid. The goal is to perform all computation within PostgreSQL for production purposes.
Is this possible?

geozelot · Answer

A naive IDW implementation, from the top of my head:

UPDATE <points> AS itp
  SET  "Z" = (
    SELECT SUM(z/d)/SUM(1/d)
    FROM   (
      SELECT smpl."Z" as z,
             ST_Distance(itp.geom, smpl.geom)^<P-value> AS d
      FROM   <point> AS smpl
      ORDER BY
             itp.geom <-> smpl.geom
      WHERE  smpl."Z" IS NOT NULL
      LIMIT  <sample_size>
    ) sq
    WHERE  ipt."Z" IS NULL
  ) q
;

Where

<P-value> is the factor applied to the inverse distance; you probably want to stay between 1.0 and 2.0
<sample_size> the amount of (known) sample points that needs to be considered in the (k) nearest neighborhood

It's been a while since I had a look at the IDW algorithm, so this may need tweaking; the subquery sq returns the list of inverse distances (when used as 1/<distance>^p), so you can go from there if you need to alter the actual IDW sums.

Spatial interpolation in PostGIS without outputting Raster

One Answer

Add your own answers!

Ask a Question