Hi,
I would like to calculate surface area (area when elevation is taken into
account) for various regions across South America.
The following code example illustrates what I am doing:
require(maptools)
require(raster)
require(rgeos)
bb <- SpatialPoints(rbind(c(-85,15),c(-31,-57)))
r <- raster(bb,nrow=100,ncol=100)
values(r) <- abs(rnorm(10000,mean=500,sd=100))
proj4string(r) <- CRS('+proj=longlat +datum=WGS84')
#Create an arbitrary polygon
poly <-
gConvexHull(SpatialPoints(rbind(c(-70,5),c(-70,-15),c(-45,-20),c(-45,-7))))
proj4string(poly) <- CRS('+proj=longlat +datum=WGS84')
#subset raster to polygon
rsub <- mask(r,poly)
#what does it look like?
data(wrld_simpl)
plot(rsub)
plot(wrld_simpl,add=T)
plot(poly,add=T,border='blue')
#Calculate surface area
z <- as.matrix(rsub)
surfaceArea(z, cellx=res(rsub)[1],celly=res(rsub)[2])
Is this the proper way to do this? My main concern is that my coordinates
are in longlat, but my elevation data (from worldclim.org) is in meters. It
seems like the correct thing to do would be to transform the elevation
raster and polygon into UTM so as to have everything in meters, but how do
you do that when the area spans across multiple UTM zones?
Any advice welcome! Thank you!
-Pascal
--
Pascal Title, MSc.
PhD student, Rabosky Lab <http://www-personal.umich.edu/~drabosky/Home.html>
Dept of Ecology and Evolutionary Biology
University of Michigan
pti...@umich.edu
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
