GIS in R

Libraries

library(raster); 
  #For raster-based loading, calculations, and mapping 
library(maptools);
  #For shapefile-based loading and mapping
library(rgdal);
  #For reprojecting pesky data

• Loading and plotting point files
• Loading and plotting rasters
• Loading and plotting shapefiles

It's as simple as loading a .csv file with a header.

shermanSquirrels <- read.csv(
  file = "./ShermanSquirrels.csv", header = T);

The format of a point file is pretty universal.

head(shermanSquirrels);

        Species Longitude Latitude
1 Sciurus niger -83.36905 31.40090
2 Sciurus niger -82.60067 28.46997
3 Sciurus niger -81.46268 28.71139
4 Sciurus niger -81.18965 28.38080
5 Sciurus niger -81.18860 28.38107
6 Sciurus niger -81.18832 28.38104

Not dissimilar from loading .csv files.

altitude <- raster(x = "./altitude.asc");

Checking out the guts of the file

altitude;

class       : RasterLayer 
dimensions  : 3600, 8640, 31104000  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -60, 90.00001  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /Users/HannahOwens/Dropbox/ENMSeminar/Labs:Homeworks/Lab2/Lab2Data/altitude.asc 
names       : altitude 

Defining coordinate reference system.

crs(altitude) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
altitude;

class       : RasterLayer 
dimensions  : 3600, 8640, 31104000  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -60, 90.00001  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/HannahOwens/Dropbox/ENMSeminar/Labs:Homeworks/Lab2/Lab2Data/altitude.asc 
names       : altitude 

plot(altitude, main="Altitude");

plot(altitude, main="Altitude", xlim = c(-88.5,-79), ylim = c(24, 32));
points(shermanSquirrels[,2:3], pch = 16, cex = 2.5)

Projections are important!

PROJCS[“NAD_1983_HARN_StatePlane_Florida_West_FIPS_0902_Feet”, GEOGCS[“GCS_North_American_1983_HARN”, DATUM[“D_North_American_1983_HARN”, SPHEROID[“GRS_1980”,6378137.0,298.257222101]], PRIMEM[“Greenwich”,0.0], UNIT[“Degree”,0.0174532925199433]], PROJECTION[“Transverse_Mercator”], PARAMETER[“False_Easting”,656166.6666666665], PARAMETER[“False_Northing”,0.0], PARAMETER[“Central_Meridian”,-82.0], PARAMETER[“Scale_Factor”,0.9999411764705882], PARAMETER[“Latitude_Of_Origin”,24.33333333333333], UNIT[“Foot_US”,0.3048006096012192]]

florida1 <- readOGR("./FloridaBound1/FLORIDA.shp");

OGR data source with driver: ESRI Shapefile 
Source: "/Users/HannahOwens/Dropbox/ENMSeminar/Labs:Homeworks/Lab2/Lab2Data/FloridaBound1/FLORIDA.shp", layer: "FLORIDA"
with 1 features
It has 4 fields
Integer64 fields read as strings:  OBJECTID 

florida1;

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -1110383, 1298633, 77020.92, 2466420  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.9999411764705882 +x_0=199999.9999999999 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs 
variables   : 4
names       : OBJECTID,  DATESTAMP,  SHAPE_AREA, SHAPE_LEN 
value       :        1, 2000/05/16, 1.56987e+12,  28344148 

plot(altitude, main="Altitude", xlim = c(-88.5,-79), ylim = c(24, 32));
plot(florida1, add = TRUE);

Oh yeah, they're at different projections. Time to warp Florida!

floridaReprojected <- spTransform(florida1, CRSobj = crs(altitude));

plot(altitude, main="Altitude", xlim = c(-88.5,-79), ylim = c(24, 32));
plot(floridaReprojected, add = TRUE);
points(shermanSquirrels[,2:3], pch = 16, cex = 2, add = TRUE);

• Masking rasters to shapefiles
• Deleting points outside of a shapefile
• Making a map

plot(altitude, main="Altitude", xlim = c(-88.5,-79), ylim = c(24, 32));
plot(floridaReprojected, add = TRUE);

altCroppedFlorida <- mask(altitude, mask = floridaReprojected)
plot(altCroppedFlorida, main="Altitude", xlim = c(-88.5,-79), ylim = c(24, 32));
plot(floridaReprojected, add = TRUE);

extent(altCroppedFlorida);

class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -60 
ymax        : 90.00001 

extent(floridaReprojected);

class       : Extent 
xmin        : -87.62601 
xmax        : -80.03095 
ymin        : 24.54522 
ymax        : 30.99702 

plot(altCroppedFlorida);

altCroppedFlorida <- crop(altCroppedFlorida, extent(floridaReprojected));
plot(altCroppedFlorida);

writeRaster(altCroppedFlorida, "croppedFlorida", format = "ascii", overwrite = TRUE);

plot(floridaReprojected, xlim = c(-88, -79), ylim = c(24, 35));
points(shermanSquirrels[,2:3], pch = 16, cex = 2);

alt <- extract(altCroppedFlorida, shermanSquirrels[2:3]);
shermanSquirrelsAlt <- cbind(shermanSquirrels, alt);
head(shermanSquirrelsAlt);

        Species Longitude Latitude alt
1 Sciurus niger -83.36905 31.40090  NA
2 Sciurus niger -82.60067 28.46997  17
3 Sciurus niger -81.46268 28.71139  24
4 Sciurus niger -81.18965 28.38080  22
5 Sciurus niger -81.18860 28.38107  22
6 Sciurus niger -81.18832 28.38104  22

nrow(shermanSquirrelsAlt);

[1] 35

shermanFlorida <- shermanSquirrelsAlt[
  complete.cases(shermanSquirrelsAlt),];
nrow(shermanFlorida);

[1] 29

plot(floridaReprojected, xlim = c(-88, -79), ylim = c(24, 35));
points(shermanSquirrels[,2:3], pch = 16, cex = 2, col = "red");
points(shermanFlorida[,2:3], pch = 16, cex = 2);

write.csv(shermanFlorida, file = "shermanSquirrelsFlorida.csv", row.names = F);

plot(altCroppedFlorida);
plot(floridaReprojected, add = TRUE);
points(shermanFlorida[,2:3], pch = 16, cex = 2);

png(filename="shermanSquirrelsFloridaMap.png");
plot(altCroppedFlorida);
plot(floridaReprojected, add = TRUE);
points(shermanFlorida[,2:3], pch = 16, cex = 2);
dev.off();

quartz_off_screen 
                2