Spatial data extraction around buffered points in R

[This article was first published on Ecology in silico, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Quantifying spatial data (e.g. land cover) around points can be done in a variety of ways, some of which require considerable amounts of patience, clicking around, and/or cash for a license.
Here’s a bit of code that I cobbled together to quickly extract land cover data from the National Land Cover Database for buffered regions around points (e.g. small ponds, point-count locations, etc.), in this case, U.S. capitals.

Here, extract_cover() is a function to do the extraction (with the help of the raster package), and extraction.R makes a parallel call to the function using the doMC package:

extract_cover.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
<span class="line"><span class="c1"># Function to extract the data in R instead of Arc</span>
</span><span class="line"><span class="c1"># inputs: year, buffer size (in meters), point data file, </span>
</span><span class="line"><span class="c1">#   cover data directory, and output file directory</span>
</span><span class="line"><span class="c1"># output: csv file with site id, cover type, and % in buffer</span>
</span><span class="line">
</span><span class="line">extract_cover <span class="o"><-</span> <span class="kr">function</span><span class="p">(</span>year<span class="p">,</span> buffer<span class="p">,</span>
</span><span class="line">                          point_d <span class="o">=</span> <span class="s">"sites.csv"</span><span class="p">,</span>
</span><span class="line">                          data_dir<span class="o">=</span><span class="s">"NLCD_data"</span><span class="p">,</span>
</span><span class="line">                          write_dir<span class="o">=</span><span class="s">"extracted"</span><span class="p">){</span>
</span><span class="line">  require<span class="p">(</span>raster<span class="p">)</span>
</span><span class="line">  require<span class="p">(</span>rgdal<span class="p">)</span>
</span><span class="line">  require<span class="p">(</span>stringr<span class="p">)</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># load cover data </span>
</span><span class="line">  filename <span class="o"><-</span> paste<span class="p">(</span>data_dir<span class="p">,</span> <span class="s">"/nlcd_"</span><span class="p">,</span>
</span><span class="line">                    year<span class="p">,</span>
</span><span class="line">                    <span class="s">"_landcover_2011_edition_2014_03_31.img"</span><span class="p">,</span>
</span><span class="line">                    sep<span class="o">=</span><span class="s">""</span><span class="p">)</span>
</span><span class="line">  NLCD <span class="o"><-</span> raster<span class="p">(</span>filename<span class="p">)</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># load site data</span>
</span><span class="line">  sites <span class="o"><-</span> read.csv<span class="p">(</span>point_d<span class="p">,</span> header<span class="o">=</span><span class="k-Variable">T</span><span class="p">)</span>
</span><span class="line">  coords <span class="o"><-</span> sites<span class="p">[,</span> c<span class="p">(</span><span class="s">"longitude"</span><span class="p">,</span> <span class="s">"latitude"</span><span class="p">)]</span>
</span><span class="line">
</span><span class="line">  <span class="c1">#convert lat/lon to appropriate projection</span>
</span><span class="line">  names<span class="p">(</span>coords<span class="p">)</span> <span class="o"><-</span> c<span class="p">(</span><span class="s">"x"</span><span class="p">,</span> <span class="s">"y"</span><span class="p">)</span>
</span><span class="line">  coordinates<span class="p">(</span>coords<span class="p">)</span> <span class="o"><-</span> <span class="o">~</span>x <span class="o">+</span> y
</span><span class="line">  proj4string<span class="p">(</span>coords<span class="p">)</span> <span class="o"><-</span> CRS<span class="p">(</span><span class="s">"+proj=longlat +ellps=WGS84 +datum=WGS84"</span><span class="p">)</span>
</span><span class="line">  crs_args <span class="o"><-</span> NLCD<span class="o">@</span>crs<span class="o">@</span>projargs
</span><span class="line">  sites_transformed <span class="o"><-</span> spTransform<span class="p">(</span>coords<span class="p">,</span> CRS<span class="p">(</span>crs_args<span class="p">))</span>
</span><span class="line">
</span><span class="line">  <span class="c1">#extract land cover data for each point, given buffer size</span>
</span><span class="line">  Landcover <span class="o"><-</span> extract<span class="p">(</span>NLCD<span class="p">,</span> sites_transformed<span class="p">,</span> buffer<span class="o">=</span>buffer<span class="p">)</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># summarize each site's data by proportion of each cover type</span>
</span><span class="line">  summ <span class="o"><-</span> lapply<span class="p">(</span>Landcover<span class="p">,</span> <span class="kr">function</span><span class="p">(</span>x<span class="p">){</span>
</span><span class="line">    prop.table<span class="p">(</span>table<span class="p">(</span>x<span class="p">))</span>
</span><span class="line">  <span class="p">}</span>
</span><span class="line">  <span class="p">)</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># generate land cover number to name conversions</span>
</span><span class="line">  num.codes <span class="o"><-</span> unique<span class="p">(</span>unlist<span class="p">(</span>Landcover<span class="p">))</span>
</span><span class="line">  cover.names <span class="o"><-</span> NLCD<span class="o">@</span>data<span class="o">@</span>attributes<span class="p">[[</span><span class="m">1</span><span class="p">]]</span><span class="o">$</span>Land.Cover.Class<span class="p">[</span>num.codes <span class="o">+</span> <span class="m">1</span><span class="p">]</span>
</span><span class="line">  levels<span class="p">(</span>cover.names<span class="p">)[</span><span class="m">1</span><span class="p">]</span> <span class="o"><-</span> <span class="kc">NA</span> <span class="c1"># first level is ""</span>
</span><span class="line">  conversions <span class="o"><-</span> data.frame<span class="p">(</span>num.codes<span class="p">,</span> cover.names<span class="p">)</span>
</span><span class="line">  conversions <span class="o"><-</span> na.omit<span class="p">(</span>conversions<span class="p">)</span>
</span><span class="line">  conversions <span class="o"><-</span> conversions<span class="p">[</span>order<span class="p">(</span>conversions<span class="o">$</span>num.codes<span class="p">),]</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># convert to data frame</span>
</span><span class="line">  mydf <span class="o"><-</span> data.frame<span class="p">(</span>id <span class="o">=</span> rep<span class="p">(</span>sites<span class="o">$</span>id<span class="p">,</span> lapply<span class="p">(</span>summ<span class="p">,</span> length<span class="p">)),</span>
</span><span class="line">                     cover <span class="o">=</span> names<span class="p">(</span>unlist<span class="p">(</span>summ<span class="p">)),</span>
</span><span class="line">                     percent <span class="o">=</span> unlist<span class="p">(</span>summ<span class="p">)</span>
</span><span class="line">  <span class="p">)</span>
</span><span class="line">
</span><span class="line">  <span class="c1"># create cover name column</span>
</span><span class="line">  mydf<span class="o">$</span>cover2 <span class="o"><-</span> mydf<span class="o">$</span>cover
</span><span class="line">  levels<span class="p">(</span>mydf<span class="o">$</span>cover2<span class="p">)</span> <span class="o"><-</span> conversions<span class="o">$</span>cover.names
</span><span class="line">
</span><span class="line">  <span class="c1"># write output</span>
</span><span class="line">  out_name <span class="o"><-</span> paste<span class="p">(</span>write_dir<span class="p">,</span> <span class="s">"/"</span><span class="p">,</span>
</span><span class="line">                    year<span class="p">,</span> <span class="s">"_cover_"</span><span class="p">,</span> buffer<span class="p">,</span>
</span><span class="line">                    <span class="s">"_m_buffer.csv"</span><span class="p">,</span> sep<span class="o">=</span><span class="s">""</span><span class="p">)</span>
</span><span class="line">  write.csv<span class="p">(</span>mydf<span class="p">,</span> out_name<span class="p">)</span>
</span><span class="line"><span class="p">}</span>
</span>

extraction.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
<span class="line"><span class="c1"># Script to actually extract the data at various buffer distances</span>
</span><span class="line"><span class="c1"># parallelized for speed (one core per year)</span>
</span><span class="line">library<span class="p">(</span>doMC<span class="p">)</span>
</span><span class="line">
</span><span class="line">years <span class="o"><-</span> c<span class="p">(</span><span class="m">2001</span><span class="p">,</span> <span class="m">2006</span><span class="p">,</span> <span class="m">2011</span><span class="p">)</span>
</span><span class="line">nyears <span class="o"><-</span> length<span class="p">(</span>years<span class="p">)</span>
</span><span class="line">registerDoMC<span class="p">(</span>nyears<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># input vector of distances (in meters)</span>
</span><span class="line">buffer_distances <span class="o"><-</span> c<span class="p">(</span><span class="m">1000</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">foreach <span class="p">(</span>i<span class="o">=</span><span class="m">1</span><span class="o">:</span>nyears<span class="p">)</span> <span class="o">%dopar%</span> <span class="p">{</span>
</span><span class="line">  <span class="kr">for</span> <span class="p">(</span>j <span class="kr">in</span> buffer_distances<span class="p">){</span>
</span><span class="line">    extract_cover<span class="p">(</span>year <span class="o">=</span> years<span class="p">[</span>i<span class="p">],</span> buffer<span class="o">=</span>j<span class="p">)</span>
</span><span class="line">  <span class="p">}</span>
</span><span class="line"><span class="p">}</span>
</span>

Resources

To leave a comment for the author, please follow the link and comment on their blog: Ecology in silico.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)