# Finding interpolated value from raster data in PostGIS

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)`.

• 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. Commented Feb 26 at 10:40
• 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. Commented Feb 27 at 9:11
• 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? Commented Feb 27 at 11:28
• @geozelot I edited my question to add a specific example of what I want. Commented Feb 27 at 13:48

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.