Extracting Raster Values from Points in R and GRASS

Saturday, July 24, 2010

A common task in GIS analysis is to extract the value of a remotely sensed environmental variable at a point location. For instance we may wish to extract the elevation of a field plot from a digital elevation model. The elevation data is a raster (i.e. grid) and the plots are a point shapefile (or a simple text file of X, Y locations). The command for doing this in ArcGIS is =ExtractValuesToPoints= available in the Spatial Analyst package. Situations may arise where ArcGIS is not the most efficient way of extracting these values. So, here, I provide a brief overview of how to extract raster values to points in R and GRASS. ** Extract Values to Points in R This is strikingly easy is R. My work usually requires more statistical sophistication than is available in ArcGIS. As a result, I have completely switched to doing the extraction in R. I known I am going to end in R eventually, and it is easier to automate than writing a long python script in ArcGIS. *** Data required For the purpose of this exercise. All the data must be have the same spatial projection. - [[file:data/gr.asc][=gr.asc=]] :: an ESRI ASCII grid. This could also be an ArcGIS binary grid if you know how to use RGDAL. That perhaps will be another post. - [[file:data/pt.zip][=pt.shp=]] :: a point shapefile. You also need the =maptools= and =sp= packages in R. *** The Code #+begin_src R -t :results output :colnames t library(maptools) # autoloads sp gr <- readAsciiGrid("data/gr.asc") pt <- readShapePoints("data/pt.shp") overlay(gr, pt) #+end_src #+results: #+begin_example coordinates gr.asc 1497 (569292, 1224170) 6094.6080 539 (567718, 1227840) 7964.5331 1023 (564565, 1225810) -293.6599 663 (562462, 1227260) -5351.7297 69 (563675, 1229710) -2923.9394 716 (563062, 1227180) -4241.8255 339 (567636, 1228780) 4781.6488 509 (561722, 1227870) -3958.7312 805 (560981, 1226690) 139.9091 155 (560884, 1229220) -2719.7353 #+end_example That is it. Fast, and easy. ** Extracting Values in GRASS Extracting raster values in GRASS is somewhat faster than in R, but it takes a little bit more planning in that you have to explicitly create the column that the raster values will go into. *** Data Required - =gr= : A GRASS grid - =pt= : A GRASS point dataset *** The Code The basic flow of this is that you create an empty column in the point dataset with the right data type (i.e. =varchar(10)= string of length 10, =double precision= floating point numbers, =int= integers). Then, fill the column with the raster values. #+begin_src sh :results: silent v.db.addcol map=pt columns="grval double precision" v.what.rast vector=pt raster=gr column=grval #+end_src