The function "extract" (from raster package) can extract raster values to points. So all you need is to convert your polygons vertices to points. The function 'as' (from package sf) can do this easily. One example:

library(raster)

library(sf)

library(rgdal)

#Read your polygons shapefile

p

Thanks a lot Tiago Gandra ! I followed your example and could extract the z-values for each point. However, I am still struggling to add the z-values back to the multipolygons, which look like this:

MULTIPOLYGON Z (((2683159 1246484 0, 2683168 12... MULTIPOLYGON Z (((2683439 1246444 0, 2683414 12... MULTIPOLYGON Z (((2683479 1246404 0, 2683452 12... MULTIPOLYGON Z (((2683159 1246484 0, 2683149 12... MULTIPOLYGON Z (((2683419 1246324 0, 2683376 12...

BTW: I have figured out that the polygon to polyline and point conversion can also be done with sf like this:

p

Michael Strupler , maybe you can add the polygon ID to each of the points and convert it to dataframe.

Then you rebuilt the polygons using the function "spatialpolygons".

Here you can find an example: https://gis.stackexchange.com/questions/171124/data-frame-to-spatialpolygonsdataframe-with-multiple-polygons

Thanks a lot, I will try that...