Interpolate raster value using Python

Question

I'm new in coding. I have a shapefile (points) and some raster files. My purpose is get the values from all raster to point (each point will get value from 2 or 3 nearest cell,the value on the value of cell and distance of the point to cell). I already built code for that and it is working but it takes so much time for running.

Could you optimize the code to make it run faster?

import rioxarray as rxr
import geopandas as gpd   

### definition extract_raster:
def extract_raster(points_gdf, raster_xarray):
    extracted_values = []
    for _, point in points_gdf.iterrows():
        value = raster_xarray.interp(x=point.geometry.x, y=point.geometry.y)
        extracted_values.append(value.values)
    
    column_name = 'raw'
    points_gdf[column_name] = extracted_values
    points_gdf[column_name] = points_gdf[column_name].astype(float)
    return points_gdf
# import data

lidar_chm_path = r"D:\CHM\CHM.tif"
lidar_chm_xr = rxr.open_rasterio(lidar_chm_path, masked=True).squeeze()

NFI_path = r"D:\SHP_file\Final_NFI.shp"
NFI_data = gpd.read_file(NFI_path)


##CHM
extract_raster(NFI_data, lidar_chm_xr)
NFI_data.rename(columns={'raw': 'CHM'}, inplace=True)

Answer 1

document the types

def extract_raster(points_gdf, raster_xarray):

Both parameters have needlessly long names. Use optional type hinting to explain to the reader that points is a geo dataframe, and raster is an xarray.

kw defaults

Consider tacking on ... , column_name: str = 'raw'): to the signature.

But now that we're promoting a local temp var to become part of the Public API, the documentation burden is greater, so maybe go with output_column_name or some similarly informative identifier.

---

comprehensions

    extracted_values = []
    for _, point in points_gdf.iterrows():
        value = raster_xarray.interp(x=point.geometry.x, y=point.geometry.y)
        extracted_values.append(value.values)

Consider using a two-line list comprehension for those five lines.

---

inplace

... .rename( ... , inplace=True)

Don't do it. Shake the habit. Here is why.

---

speed

This submission is about performance, but it did not include any profiling data, characterize the inputs of the target workload, nor mention any timing measurements. It only lamented that "it takes so much time for running."

optimize the code

There just isn't a lot of code there, so there isn't a lot of flexibility for choosing slightly different approaches. With some profiling data in hand we might delve into options for open_rasterio() or read_file() to see if we could ask the more expensive one to do less work.

If we had a reprex we might talk about how to let numpy or interp() do the looping over points, but we don't see one in the OP.

Converting a loop to a list comprehension will improve source code readability, but won't substantially affect the elapsed running time.

It's unclear how the extracted rasters deliver business value. If similar value could be delivered with subsampled data, say a 4x or 16x reduction in the number of returned points, then go ahead and process less data in less time. We're interpolating anyway, so maybe focusing attention on local maxima and minima would go a long way toward reducing amount of data without much sacrifice of accuracy. Curve fitting, perhaps with splines, may allow significant compression.