Extracting Raster Values from Points in R and GRASS

Monday, July 26, 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.

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.

pt.shp

a point shapefile.

You also need the maptools and sp packages.

The Code

library(maptools) # autoloads sp gr <- readAsciiGrid("www.unc.edu/~toddjobe/data/gr.asc") pt <- readShapePoints("www.unc.edu/~toddjobe/data/pt.shp") overlay(gr, pt)

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

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.

v.db.addcol map=pt columns="grval double precision" v.what.rast vector=pt raster=gr column=grval

Tmp

I use org-mode in Emacs for my work logs, for writing manuscripts, for taking notes, and, increasingly, for writing my code. I also use org-mode to write this blog. Today, I'll share how I publish a blog like mine hosted on Google's blogger.com using org-mode . Users unfamiliar with emacs can get a quick tour here. Emacs users unfamiliar with org-mode will be convinced to try it by watching this video. Org-mode increased my work efficiency greatly.

I keep my blog in a single org file: blog.org. Here is a snippet:

Tmp

I use org-mode in Emacs for my work logs, for writing manuscripts, for taking notes, and, increasingly, for writing my code. I also use org-mode to write this blog. Today, I'll share how I publish a blog like mine hosted on Google's blogger.com using org-mode . Users unfamiliar with emacs can get a quick tour here. Emacs users unfamiliar with org-mode will be convinced to try it by watching this video. Org-mode increased my work efficiency greatly.

I keep my blog in a single org file: blog.org. Here is a snippet:

#+TITLE:     blog.org
#+AUTHOR:    R. Todd Jobe
#+EMAIL:     toddjobe@unc.edu
#+DATE:      %D
#+DESCRIPTION: See Todd Jobe's thoughts about computers and ecology. R code. Arc-GIS and GRASS code. Bash scripts. Matlab code.
#+OPTIONS:   H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
#+OPTIONS:   TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc
#+INFOJS_OPT: view:nil toc:nil ltoc:t mouse:underline buttons:0 path:http://orgmode.org/org-info.js
#+EXPORT_SELECT_TAGS: export
#+EXPORT_EXCLUDE_TAGS: noexport
#+TODO: I(i) S(s) R(r!) | Q(q!) P(p!/!) C(@/!) (todo)

* Q Extracting Raster Values from Points in R and GRASS
  :PROPERTIES:
  :Tags: R, GRASS, GIS
  :END:
  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
** 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 -t -w 66 :results: silent
   v.db.addcol map=pt columns="grval double precision"
   v.what.rast vector=pt raster=gr column=grval
#+end_src

* S Posting to blogger using org-mode
   I use [[http://orgmode.org/][org-mode]] in Emacs for my work logs, for writing manuscripts,
   for taking notes, and, increasingly, for writing my code.  I also
   use org-mode to write this blog.  Today, I'll share how I publish a
   blog like mine hosted on Google's [[http://blogger.com/][blogger.com]] using org-mode .
   Users unfamiliar with emacs can get a quick tour [[http://www.gnu.org/software/emacs/tour/][here]].  Emacs users
   unfamiliar with org-mode will be convinced to try it by watching
   [[http://www.youtube.com/watch?v%3DoJTwQvgfgMM][this]] video.  Org-mode increased my work efficiency greatly.
   
   I keep my blog in a single org file: ~blog.org~.  Here is a
   snippet: 

Each blog entry
   is a TODO entry.  As I 
   blog
   ~blog.org~. 
   My [[http://toddjobe.ecology@blogspot.com][blog]] is hosted on [[http://blogger.com][blogger.com]].  This service used to allow me to
   FTP my blog directly from personal website:
   [[http://www.unc.edu/~toddjobe][www.unc.edu/~toddjobe/]].  Google recently ceased support for FTPing
   since less than 0.5% of their users employed this feature.  So, now
   I have to upload directly to 
*** Export region as html
*** Use Google CL and the clipboard
* S Changing the font in Carbon Emacs

  My current monospaced font of choice is [[http://dejavu-fonts.org/wiki/index.php%3Ftitle%3DMain_Page][Deja Vu Sans Mono]].  It is
  pretty easy to change the default font in Carbon Emacs.  It is not
  so easy to change the default font on Windows- or Linux-based
  Emacs22, which is the current release version.  It will be
  relatively straightforward in Emacs23 via the ~xft~ package.  But,
  for your Mac folks out there...

** Installing a Font on Mac OS X
   1. Download the font (Mac will accept almost any format).
   2. Drag the font file into ~/System/Fonts~

** Figure out the long name of your font face in Carbon Emacs
   - M-x ~mac-font-panel-mode~.  Then pick your font
   - M-x ~describe-font~ and copy the long name.

** Set the default font
   There are a couple of ways to do this.
*** Customization
    - M-x ~customize-face~ RET ~default~ RET
    - Fill in the fields of the font with your 
*** ~.emacs.d~ file
#+BEGIN_SRC elisp 
(set-default-font "-apple-inconsolata-medium-r-normal--11-110-72-72-m-110-iso10646-1")
#+END_SRC

** Monospaced fonts for the uninitiated
  Programmers tend to use different fonts for code than what we see
  in books an on the web.  Word processors typically use proportional
  fonts, where the horizontal space allowed for narrow letters like
  "l" is less than the space allowed for wide letters like "W".  In
  code, we want the columns of each line to be aligned.  So,
  programmers tend to use monospaced fonts, where each letter takes
  up the same amount of horizontal space.  The most commonly known
  fixed width font is Courier.  Another common fixed-width font for
  Windows and Microsoft products in Consolas.  [[http://www.codeproject.com/KB/work/FontSurvey.aspx][The Code Project]] has
  an exhaustive list of monospaced fonts.

  One of the key attributes to look for in a monospaced font is the
  distinction between the number one and the lowercase letter "l" and
  the number zero and the uppercase "O". When you are trying to read
  variable names that have mixtures of numbers and letters, you will be
  thankful for a font that distinguishes between these.
* I Spelling and Grammar checking in Emacs
  - Grammar checking is pretty easy with ~ispell~
  - The ~diction~ command line utility is a Unix classic
  - An minor-mode for calling ~diction~ from Emacs is [[http://ftp2.de.freebsd.org/pub/emacs/emacs-lisp/archive/diction.el][~diction~]]
  - Use ~diction-buffer~ or ~diction-region~ commands
* I Converting R tables to data frames.
  as.matrix doesn't work
  class(tbl) <- c("matrix")

Each blog entry
   is a TODO entry.  As I 
   blog
   ~blog.org~. 
   My [[http://toddjobe.ecology@blogspot.com][blog]] is hosted on [[http://blogger.com][blogger.com]].  This service used to allow me to
   FTP my blog directly from personal website:
   [[http://www.unc.edu/~toddjobe][www.unc.edu/~toddjobe/]].  Google recently ceased support for FTPing
   since less than 0.5% of their users employed this feature.  So, now
   I have to upload directly to 
*** Export region as html
*** Use Google CL and the clipboard
* I Spelling and Grammar checking in Emacs
  - Grammar checking is pretty easy with ~ispell~
  - The ~diction~ command line utility is a Unix classic
  - An minor-mode for calling ~diction~ from Emacs is [[http://ftp2.de.freebsd.org/pub/emacs/emacs-lisp/archive/diction.el][~diction~]]
  - Use ~diction-buffer~ or ~diction-region~ commands

The first thing to notice is that this org-file has custom TODO states

Tmp

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.

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.

pt.shp

a point shapefile.

You also need the maptools and sp packages in R.

The Code

library(maptools) # autoloads sp gr <- readAsciiGrid("www.unc.edu/~toddjobe/data/gr.asc") pt <- readShapePoints("www.unc.edu/~toddjobe/data/pt.shp") overlay(gr, pt)

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

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.

v.db.addcol map=pt columns="grval double precision" v.what.rast vector=pt raster=gr column=grval

Tmp

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.

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.

pt.shp

a point shapefile.

You also need the maptools and sp packages in R.

The Code

library(maptools) # autoloads sp gr <- readAsciiGrid("www.unc.edu/~toddjobe/data/gr.asc") pt <- readShapePoints("www.unc.edu/~toddjobe/data/pt.shp") overlay(gr, pt)

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

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.

v.db.addcol map=pt columns="grval double precision"
v.what.rast vector=pt raster=gr column=grval

tmp

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.

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.

pt.shp

a point shapefile.

You also need the maptools and sp packages in R.

The Code

library(maptools) # autoloads sp gr <- readAsciiGrid("www.unc.edu/~toddjobe/data/gr.asc") pt <- readShapePoints("www.unc.edu/~toddjobe/data/pt.shp") overlay(gr, pt)

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

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.

v.db.addcol map=pt columns="grval double precision"
v.what.rast vector=pt raster=gr column=grval

New post

Tmp