Load libraries and grab tree cover data: https://catalog.data.gov/dataset/nlcd-2011-percent-tree-canopy-cartographic
Tree canopy cover ranges from 0-100%.
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.3.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.3.2
## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4
library(raster)
x<-raster("nlcd_2011_USFS_tree_canopy_2011_edition_2014_03_31/cartographic_product/nlcd2011_usfs_treecanopy_cartographic_3-31-2014.img")
plot(x)
Ok, now load in dataset. You just need a dataset with lat, lon.
library(data.table)
## Warning: package 'data.table' was built under R version 3.3.2
##
## Attaching package: 'data.table'
## The following object is masked from 'package:raster':
##
## shift
b<-fread("https://raw.githubusercontent.com/adnguyen/adnguyen.github.io/master/demos/20160517_ANBE_ant_sampling.csv",skip=6) #grab data from online
#b<-read.csv("20160517_ANBE_ant_sampling.csv",skip=6) # if data is local, you can use this
b<-b[-87,] # getting rid of row that does not have lat lon
coords <- b[, c("lon", "lat")] # grabbing just the coordinates
names(coords) <- c("x", "y")
coordinates(coords) <- ~x + y
proj4string(coords) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
crs_args <- x@crs@projargs
sites_transformed <- spTransform(coords, CRS(crs_args))
#extract land cover data for each point, given buffer size
#Landcover <- extract(x, sites_transformed, buffer=30)
Landcover <- extract(x, sites_transformed, buffer=20)
Landcover[1:10] # you get a list
## [[1]]
## value
## 74
##
## [[2]]
## value
## 74
##
## [[3]]
## value
## 74
##
## [[4]]
## value
## 23
##
## [[5]]
## value
## 23
##
## [[6]]
## value
## 23
##
## [[7]]
## value
## 0
##
## [[8]]
## value
## 0
##
## [[9]]
## value
## 0
##
## [[10]]
## [1] 0 0
b$tree_canopy_20<-unlist(lapply(Landcover, `[[`, 1))
### just grabbing a smaller buffer radius
Landcover2 <- extract(x, sites_transformed, buffer=.25)
Landcover2[1:10]
## [[1]]
## [1] 74
##
## [[2]]
## [1] 74
##
## [[3]]
## [1] 74
##
## [[4]]
## [1] 23
##
## [[5]]
## [1] 23
##
## [[6]]
## [1] 23
##
## [[7]]
## [1] 0
##
## [[8]]
## [1] 0
##
## [[9]]
## [1] 0
##
## [[10]]
## [1] 0
length(unlist(Landcover2)) # make sure this number is the same length as your data frame
## [1] 114
b$tree_canopy<-unlist(Landcover2) # link up tree cover to original data
#names(b)
#b
#write.csv(b,"test.csv") # write out dataset
Depending on buffer radius, you’ll get a few values per site. So I just took the first value.
sessionInfo()
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X 10.12.2 (Sierra)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] data.table_1.10.0 raster_2.5-8 rgdal_1.2-5 sp_1.2-4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.9 lattice_0.20-34 digest_0.6.12 rprojroot_1.2
## [5] grid_3.3.1 backports_1.0.5 magrittr_1.5 evaluate_0.10
## [9] stringi_1.1.2 curl_2.3 rmarkdown_1.3 tools_3.3.1
## [13] stringr_1.2.0 yaml_2.1.14 htmltools_0.3.5 knitr_1.15.1