Extracting raster values at points using Open Source GIS?

QGIS "Point Sampling Tool" should be the plugin you're looking for.

Here's a detailed description of how to use it: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Update based on Paolo's comment:

the plugin is not the only solution, and not always the easiest solution anymore. An alternative solution is the Saga function ‘Add raster values to point’ in the processing toolbox. See for details http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

I was having problems with the QGIS and SAGA GUI tools mentioned in this thread (Raster values to points was failing for some reason and throwing unhelpful errors and the GRASS v.sample created a whole new layer which was not helpful). After failing with the GUI tools for a while, I tried doing this in the Field Calculator. It worked quite well and I was able to control the process a little better than the GUIs allow, and make some other calculations along the way.

Say you have a layer named pts and another named rast, both in the same coordinate system. You'd like to sample rast at each X,Y pair represented in pts.

If you haven't used the Field Calculator before, it's pretty simple. You will enter your calculation in the "Expression" box, and Q gives you a number of variables and operations in the middle column, with help text on in the right column. I'll break this process into four steps:

1. Open the attribute table of the pts layer you'd like to sample with.

2. Once you are in the Field Calculator dialog, choose whether you'd like to Create a new field or Modify an existing field in your pts layer.

3. Next, build an expression to fill the new or existing pts attribute column. You might start by modifying the expression code that worked for me:

raster_value('rast', 1, make_point($x, $y))

4. You must supply raster_value() with a raster layer name 'rast', a band number 1, and the point geometry at make_point(). $x and $y are geometry variables dependent on the location of the point in each row of the attribute table.

This method also allows arithmetic operations like subtracting the value of another raster layer called other_rast from rast, which saved me a bunch of time over the GUI tools. Example below:

raster_value('rast', 1, make_point($x, $y)) - raster_value('other_rast', 1, make_point($x, $y))

Note again that the three layers pts, rast, and other_rast must be in the same coordinate system for this method to work.

In PostGIS 2.0 you can do:

SELECT ST_Value(rast, geom) val
FROM yourrastertabe, yourpointtable
WHERE ST_Intersects(rast, geom)

Make sure your raster is tiled very small when you load it (-t 10x10 with the loader).