Finding interpolated value from raster data in PostGIS

Question

I have a series of values that were computed for each cell on a grid on a given territory.

Those values has been stored in a raster field in PostGIS (there are no missing data).

I want to run a query that will give me the value for a given coordinates.

Now, the grid can be fairly large, and I want to get an interpolated value, depending on where does the coordinates fall inside the cell.

My problem is that I can only ask PostgreSQL for the nearest values, or for the nearest neighborhood.

SELECT ST_Neighborhood(tile.rast, point, 1, 1)
FROM geodata_catchmentareatile as tile
CROSS JOIN ST_Transform(ST_Point(-1.45654, 46.42324, 4326), 2154) AS point
WHERE ST_Intersects(tile.rast, point);

=> {{6895,10917,18151},{5292,13371,14158},{6895,10939,11704}}

Here I get a grid of values for the pixels around my coordinates, but if I move said coordinates from just a few meters in any direction, I might fall in another cell and then get a huge jump of value, whereas I would like something more continuous.

I don't have the slightest clue as to how to achieve this. Any pointers?

-- edit --

To be more specific, let's imagine the raster field stored in PostGIS contains the following data...

[
 (100, 100), 1000,
 (100, 110), 1500,
 (110, 100), 4000,
 (110, 110), 250,
]

...where (100, 100) are the actual coordinates for a data point (and where the coordinates system is in meters), and 1000 is the measurement at this point. I want to be able to query an interpolated value at a specific position, e.g (105, 108).

• RoperMaps: Would taking the average of the non-null values returned not achieve this? If you want greater smoothing increase distanceX and distanceY. This sometimes called a moving or sliding average. I think soem of the confusion might be about the difference between interpolation and smoothing. I think you want smoothing.
• Thibault J: Thank you for your comment. I don't think what I want can be achieved using a rolling average or smoothing, since it still would not take into account where exactly is my query point located relative to the data grid.
• geozelot: Can you define "continuous" in this context? A new cell has new values, or a new set of neighboring values - that's the point of it. Do you want those values weighted for the actual distances (e.g. of their centroids) to those inout coordinates?
• Thibault J: @geozelot I edited my question to add a specific example of what I want.

Answer 1

Ok, I think I finally managed to find a solution.

I start by creating a squared buffer around my coordinates, using ST_Envelope and ST_Buffer. Then I clip the raster with ST_Clip and then I extract a table of (x, y, val) using ST_PixelAsPoints.

Here is the python code:

from django.contrib.gis.geos import Point
# Raster data is stored in Lambert93, and I want to query a Lat / Lng coordinates
# The grid size is 20m
EPSG_WGS84 = 4326
EPSG_LAMB93 = 2154
lng_lat = Point(float(lng), float(lat), srid=EPSG_WGS84)
lamb93_coords = lng_lat.transform(EPSG_LAMB93, clone=True)
query = """
SELECT ST_X(geom), ST_Y(geom), val
FROM (
 SELECT
 (ST_PixelAsPoints(
 ST_Clip(
 tiles.rast,
 envelope
 )
 )).*
 FROM
 geodata_catchmentareatile AS tiles
 CROSS JOIN
 ST_Transform(ST_Point(%s, %s, 4326), 2154) AS point
 CROSS JOIN
 ST_Envelope(ST_Buffer(point, 30)) AS envelope
 WHERE
 ST_Intersects(tiles.rast, envelope)) foo;
"""
cursor.execute(query, [lng, lat])
pixels = cursor.fetchall()
coords = [(x, y) for x, y, v in pixels]
values = [v for x, y, v in pixels]
interpolated_value = griddata(
 coords, values, lamb93_coords, method="cubic"
)

--edit--

There was an issue with the previous version of this query, where there would miss data if the query point was at the edge of several tiles.