Extracting Data from Rasters

This document shows you how to extract data from rasters.

 

Getting The Libraries

First, I’ll load in some packages to get the ability to work with raster data and to load in the Arapatus attenuatus data set (it is part of the default gstudio package).

require(raster)
## Loading required package: raster
## Warning: package 'raster' was built under R version 3.6.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.6.2
require(gstudio)
## Loading required package: gstudio
## Registered S3 method overwritten by 'spdep':
##   method   from
##   plot.mst ape
## Registered S3 method overwritten by 'gstudio':
##   method               from    
##   as.data.frame.genind adegenet

Loading and Cropping Rasters

We can load in the raster, and then crop it to just the are we need. These rasters were downloaded from [http://www.worldclim.org] and are much larger than the study area. This just makes it easier on the computer to not have to deal with such large areas. After cropping it, we will load in the annual precip and temperature data as well.

tmp <- raster("alt_22.tif")
e <- extent(c(-115, -109, 22, 30))
baja_alt <- crop(tmp, e)
baja_temp <- crop(raster("bio1_22.tif"), e)
baja_prec <- crop(raster("bio12_22.tif"), e)

Getting Example Data from Araptus attenuatus

Now, lets grab the Araptus data and look at the data and plot out the locations.

library(ggplot2)
data(arapat)
summary(arapat)
##       Species      Cluster      Population        ID         Latitude    
##  Cape     : 75   CBP-C :150   32     : 19   101_10A:  1   Min.   :23.08  
##  Mainland : 36   NBP-C : 84   75     : 11   101_1A :  1   1st Qu.:24.59  
##  Peninsula:252   SBP-C : 18   Const  : 11   101_2A :  1   Median :26.25  
##                  SCBP-A: 75   12     : 10   101_3A :  1   Mean   :26.25  
##                  SON-B : 36   153    : 10   101_4A :  1   3rd Qu.:27.53  
##                               157    : 10   101_5A :  1   Max.   :29.33  
##                               (Other):292   (Other):357                  
##    Longitude          LTRS          WNT            EN           EF     
##  Min.   :-114.3   01:01 :147   03:03  :108   01:01  :225   01:01 :219  
##  1st Qu.:-113.0   01:02 : 86   01:01  : 82   01:02  : 52   01:02 : 52  
##  Median :-111.5   02:02 :130   01:03  : 77   02:02  : 38   02:02 : 90  
##  Mean   :-111.7                02:02  : 62   03:03  : 22   NA's  :  2  
##  3rd Qu.:-110.5                03:04  :  8   01:03  :  7               
##  Max.   :-109.1                (Other): 15   (Other): 16               
##                                NA's   : 11   NA's   :  3               
##      ZMP           AML           ATPS          MP20    
##  01:01 : 46   08:08  : 51   05:05  :155   05:07  : 64  
##  01:02 : 51   07:07  : 42   03:03  : 69   07:07  : 53  
##  02:02 :233   07:08  : 42   09:09  : 66   18:18  : 52  
##  NA's  : 33   04:04  : 41   02:02  : 30   05:05  : 48  
##               07:09  : 22   07:09  : 14   05:06  : 22  
##               (Other):142   08:08  :  9   (Other):119  
##               NA's   : 23   (Other): 20   NA's   :  5
coords <- strata_coordinates(arapat)
##
## This no longer works UNLESS you have a google API key
## 
library(ggmap)
map <- population_map(coords, map.source = "stamen")
ggmap( map ) + geom_point(aes(x=Longitude,y=Latitude), data=coords)

Extracting Point Data

To elevation, temperature and precipitation from the rasters for each sampling location, we need to translate them into points first. I’ll first grab the coordinate data as a data.frame.

Then we can grab them using the normal functions in the sp library.

pts <- SpatialPoints(coords[, 2:3])
coords$elev <- extract(baja_alt, pts)
coords$temp <- extract(baja_temp, pts)
coords$prec <- extract(baja_prec, pts)
coords
##     Stratum Longitude Latitude elev temp prec
## 1        88 -114.2935 29.32541  681  178  143
## 11        9 -113.9449 29.01457  361  195  148
## 20       84 -113.6679 28.96651  368  197  124
## 29      175 -113.4897 28.72796  240  205  106
## 36      177 -113.9914 28.66056  177  203  120
## 46      173 -112.8698 28.40846   26  223  102
## 56      171 -113.1826 28.22308  522  195  145
## 66       89 -113.3999 28.03661  290  203  117
## 76      159 -113.3161 27.52944   87  211  123
## 85      SFr -112.9640 27.36320  305  205  107
## 94      160 -112.5296 27.40498  371  212  124
## 104     162 -112.4080 27.20280  248  220  119
## 114      12 -112.6655 27.18232  488  202  130
## 124     161 -112.9860 27.03670   29  216   91
## 134      93 -112.0461 26.94589   66  230   81
## 144     165 -112.4095 26.24905   NA   NA   NA
## 154     169 -111.3783 26.20876    6  238  134
## 164      58 -111.3547 26.01550   12  240  129
## 173     166 -112.0806 25.91409   70  224   84
## 183      64 -111.3264 25.60521  397  214  245
## 188     168 -111.2156 25.55757  354  217  235
## 198      51 -111.6006 25.34819   80  222  140
## 205   Const -111.6750 25.02470   50  217  127
## 216      77 -110.6917 24.87611   52  231  193
## 226     164 -111.5441 24.74642   35  213   95
## 236      75 -110.7460 24.58843   21  231  166
## 247     163 -110.9510 24.21150  196  219  127
## 257    ESan -110.3686 24.45879    9  237  175
## 265     153 -110.4624 24.13389   31  235  165
## 275      48 -110.2725 24.21441  117  233  217
## 285     156 -109.9890 24.04380    0  237  196
## 291     157 -110.0960 24.01950  609  208  443
## 301      73 -109.8507 24.00789  103  233  227
## 311     Aqu -110.1043 23.28550  142  226  227
## 319     Mat -110.1091 23.08984    5  234  159
## 324      98 -109.6487 23.07570   75  237  233
## 328     101 -110.5744 27.90509    8  249  244
## 337      32 -109.3270 26.63783   18  244  337
## 356     102 -109.1263 26.38014   10  245  346

Plotting Trend lines.

Cool, lets sort this by latitude

coords <- coords[order(coords$Latitude), ]

and then plot out some values to look at what is going on.

ggplot(coords, aes(x = Latitude, y = elev)) + geom_line(color = "gray") + theme_bw() + geom_text(aes(y = elev + 10, label = Stratum), color = "blue")
## Warning: Removed 1 rows containing missing values (geom_text).

ggplot(coords, aes(x = Latitude, y = temp)) + geom_line(color = "gray") + theme_bw() + geom_text(aes(y = temp + 5, label = Stratum), color = "blue")
## Warning: Removed 1 rows containing missing values (geom_text).

ggplot(coords, aes(x = Latitude, y = prec)) + geom_line(color = "gray") + theme_bw() + geom_text(aes(y = prec + 10, label = Stratum), color = "blue")
## Warning: Removed 1 rows containing missing values (geom_text).

Father, Husband, Brewer, Professor

Middle aged guy trying to keep it all together and figure out how to best navigate the world as it is. Technology geek, practitioner of fermentation sciences, researcher, biologist.

Related