Reading raster values in points gives back a generator object instead of actual value

Question

Given a raster TIF file and a point feature class in a geodatabase, I am trying to extract the values of the raster on the points and add them as a new column on the point feature class. I am using the following piece of code:

import ....
gdb_path = r'.......gdb'
pnts = gpd.read_file(gdb_path, driver='FileGDB', layer='...')
raster_path = r'.........tif'

#building list of pairs of coordinates
temp = [pnts['geometry'].x, pnts['geometry'].y]
coords = list(map(list, zip(*temp)))

with rasterio.open(raster_path) as src:
    pnts['values'] = src.sample(coords)

The code above creates a new column on the points feature class where the value for each point is a generator object sample_gen at 0x000..... instead of the actual raster value.

I could loop point by point instead, like this:

raster = rasterio.open(raster_path)
for id, row in pnts.iterrows():
    for val in raster.sample([(row['geometry'].x, row['geometry'].y)]):
        pnts.at[id, 'value'] = val[0]

This seems to work ok, but it takes a really long time (I have 4 millions of points). So, this solution is not valid for my purposes.

Could anyone shed some light on how can I make the first piece of code produce what I need? Or any other approach?

user2856 · Accepted Answer · 2019-12-19T03:11:30.900

The samples are lazily evaluated and returned as a generator so it can scale - i.e you can sample large numbers of points from a massive raster without having to load the raster, all the points and sampled values in memory.

Each element the generator returns is a list of values, one for each band in your raster. This is true even if there's only one band and thus a single sampled value for each point.

So you can't just use pnts['values'] = list(src.sample(coords)) otherwise for a single band raster you end up with something like:

print(pnts['values'])
0       [0.5706134]
1      [0.39265066]
2       [0.7231091]
3       [0.8030734]
4       [0.4549229]
           ...

So you need to grab the first element if you have a single band raster, something like:

with rasterio.open(raster_path) as src:
    pnts['values'] = [sample[0] for sample in src.sample(coords)]

print(pnts['values'])
0      0.570613
1      0.392651
2      0.723109
3      0.803073
4      0.454923
         ...

This seems to work really well (fast enough). And thanks for the clarification about the generator. — Pitrako Junior, Dec 19 '19 at 16:57

score 1 · Answer 2 · answered Dec 18 '19 at 15:10

I answered a somewhat similar question a while ago and I think it can be adapted for your example. Hopefully it runs faster than your actual solution. Also, I haven't worked with rasterio but as it is build on top of gdal, I assume you can just use gdal as well.

Here is a function that will get the indices for the points (coordinates) you want (note that the spatial reference of both the points and the raster should be the same):

import numpy


def get_indices(x, y, ox, oy, pw, ph):
    """
    Gets the row (i) and column (j) indices in an array for a given set of coordinates.
    Based on https://gis.stackexchange.com/a/92015/86131

    :param x:   array of x coordinates (longitude)
    :param y:   array of y coordinates (latitude)
    :param ox:  raster x origin
    :param oy:  raster y origin
    :param pw:  raster pixel width
    :param ph:  raster pixel height
    :return:    row (i) and column (j) indices
    """

    i = np.floor((oy-y) / ph).astype('int')
    j = np.floor((x-ox) / pw).astype('int')

    return i, j

Having this, you can open the raster, get its origin and dimensions and then pass them along with the point's x and y coordinates to the function:

import gdal

ds = gdal.Open(raster_path, 0)
xmin, xres, xskew, ymax, yskew, yres = ds.GetGeoTransform()

idx = get_indices(pnts['geometry'].x.values, pnts['geometry'].y.values, xmin, ymax, xres, -yres)

Finally, you just have to read the raster as a numpy array and index it to get the values. Note that the raster should have only a single band or you should just read one of the bands. Otherwise, you would need to get the indices for the third dimension as well.

arr = ds.ReadAsArray()
pnts['values'] = arr[idx]
del ds

"I haven't worked with rasterio" - give it a try, you won't regret it :) — user2856, Dec 19 '19 at 03:13
Thx Marcelo for your answer. I believe your suggestion would also require looping point by point, correct? (which is what I am trying to avoid, because I think that's what makes it so slow) — Pitrako Junior, Dec 19 '19 at 16:50
This leverages numpy vectorization so you are not really doing normal Python for loops. You can compare this to the other answer and your current solution and see which one runs faster. — Marcelo Villa, Dec 20 '19 at 13:15

Reading raster values in points gives back a generator object instead of actual value

2 Answers2