Week 3 — Introduction to spatial data in R
\[Ziru~Mo\] \[912922344\]
Upload this Rmarkdown file (and a word or pdf compilation) after adding R code and some text such that it does/answer the following questions.
Consider answering by writing sentences that use inline R commands like this 5 + 5 = 10.
1) read the shapefile “galap.shp” (in galap.zip) into R to create a SpatialPolygonsDataFrame.
library(sp)
library(raster)
library(rgdal)## rgdal: version: 1.3-6, (SVN revision 773)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1unzip("galap.zip")
g <- shapefile("galap.shp")2) use R code to show the class of object g
class(g)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"3) how many polygons are there in this object?
polygons(g)## class : SpatialPolygons
## features : 30
## extent : 610242.7, 918583.1, -156287.7, 185890.3 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0There are 30 polygons.
4) how many variables are there in this object; and what are their names?
ncol(g)## [1] 12There are 12 variables in this object.
names(g)## [1] "NAME_0" "NAME_1" "NAME_2" "Island" "ID"
## [6] "species" "native" "area" "elevation" "distNear"
## [11] "distSC" "areaAdj"5) make a plot of the islands (their outlines)
plot(g, col = "lightgreen", main = "Plot of Islands")
6) make a scatter plot of the number of species as a function of the size of the island
plot(g$species, g$area,
xlab = "Number of Species",
ylab = "Size of Island",
main = "Number of Species vs. Size of Island")
7) select the largest island (with code), print its name, and save the data to a shapefile
largest_island <- g[g$area == max(g$area), "Island"]
print(largest_island)## class : SpatialPolygonsDataFrame
## features : 1
## extent : 656372.2, 746962.5, -116284.8, 18536.94 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 1
## names : Island
## value : Isabelaplot(largest_island,col = "lightgreen", main = "Isabela Island")
###Save the file to shapeshile
shapefile(largest_island, filename = 'Isabela Island',overwrite=TRUE)
largest_island## class : SpatialPolygonsDataFrame
## features : 1
## extent : 656372.2, 746962.5, -116284.8, 18536.94 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 1
## names : Island
## value : Isabela8) download elevation data for Ecuador (use the “getData” function in the raster package)
Ecuador = getData('alt', country = "Ecuador", download = TRUE)
Ecuador## class : RasterLayer
## dimensions : 804, 2004, 1611216 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -91.8, -75.1, -5.1, 1.6 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : /Users/mzr/Documents/Study/ABT 181N/ECU_msk_alt.grd
## names : ECU_msk_alt
## values : -13, 6074 (min, max)9) use the crop function in the raster package to cut out the Galapagos Islands from the rest of Ecuador. Then map the elevation data of the Galapagos islands (add the island outlines). Note that the island outlines have a different CRS than the elevation data.
Showing the Initial Coordinate System
g@proj4string## CRS arguments:
## +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0Ecuador@crs## CRS arguments:
## +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0Transfer
g2 = spTransform(g, CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
g2@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0Crop The Galapagos Islands from the rest of Ecuador
gala = crop(x = Ecuador, y = g2)
gala## class : RasterLayer
## dimensions : 361, 307, 110827 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -91.8, -89.24167, -1.408333, 1.6 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : /private/var/folders/kq/147vdycn2wv68289ps0x8kbh0000gn/T/RtmpWsgSvr/raster/r_tmp_2019-01-23_201201_6424_16209.grd
## names : ECU_msk_alt
## values : -13, 1634 (min, max)Plot The Original and Cropped Shapefile
Original
plot(Ecuador, main = "Ecuador")
Cropped
plot(gala, main = "Galapagos Islands")
lines(g2)
10) compute the maximum elevation for each island. Sort the results, and show the names and average elevation of the five islands with the highest maximum elevation
elevation_max <- extract(gala, g2, fun=max, na.rm=TRUE)## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Inf
## Warning in fun(res[[i]], na.rm = na.rm): no non-missing arguments to max;
## returning -Infelevation_mean<- extract(gala, g2, fun=mean, na.rm=TRUE)
fiveisland <- data.frame(g2$Island, elevation_max, elevation_mean)
sortisland <- fiveisland[order(fiveisland$elevation_max, decreasing = TRUE),]
sortisland## g2.Island elevation_max elevation_mean
## 2 Isabela 1634 312.35675
## 3 Fernandina 1428 328.29437
## 20 San Salvador 868 189.03096
## 16 Santa Cruz 789 176.26222
## 10 San Cristóbal 643 148.79445
## 6 Pinta 556 193.28333
## 29 Santa Maria 406 170.19786
## 5 Marchena 262 109.17323
## 15 Santa Fé 190 120.50000
## 8 Española 173 75.83077
## 21 Baltra 58 30.73333
## 19 Rabida 57 57.00000
## 23 Bartholomé 21 21.00000
## 26 Eden 12 12.00000
## 1 Tortuga -Inf NaN
## 4 Las Plazas -Inf NaN
## 7 Genovesa -Inf NaN
## 9 Gardner -Inf NaN
## 11 Gardner -Inf NaN
## 12 Caldwell -Inf NaN
## 13 Enderby -Inf NaN
## 14 Campión -Inf NaN
## 17 Coamaño -Inf NaN
## 18 Pinzón -Inf NaN
## 22 Seymour -Inf NaN
## 24 Daphne Major -Inf NaN
## 25 Daphne Minor -Inf NaN
## 30 Onslow -Inf NaN
## 27 Darwin NA NA
## 28 Wolf NA NATop 5 island’s elevation
top5 <- head(sortisland,5)
top5## g2.Island elevation_max elevation_mean
## 2 Isabela 1634 312.3567
## 3 Fernandina 1428 328.2944
## 20 San Salvador 868 189.0310
## 16 Santa Cruz 789 176.2622
## 10 San Cristóbal 643 148.794511) why does the above code give so many warnings (not errors!)?
When we crop the the islands, some of the numbers are too small to show the maximum elevation which it has been counted as negative infinity. Therefore the code gives us warnings.