[R-sig-Geo] calculate surfaceArea across a broad extent in R

Michael Sumner mdsumner at gmail.com
Sun Dec 22 22:38:41 CET 2013


On Mon, Dec 23, 2013 at 7:25 AM, Pascal Title <ptitle at umich.edu> wrote:
> 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?
>

Just don't use UTM!  There are many other good choices that don't have
this zone limitation. (Folks who peddle/refer-to UTM as if it's
somehow the *only* planar projection out there have a lot to answer
for . . .)

It really depends on your study region and needs, but you probably
cannot go wrong with Lambert Azimuthal Equal Area, build your own such
as:

"+proj=laea +lon_0=-60 +lat_0=-21 +datum=WGS84"

Please explore this with your own parameter values and seek further
expert advice for your study.  There are many many choices for an
appropriate projection and many guides already exist.

Cheers, Mike.


> 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
> ptitle at umich.edu
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com



More information about the R-sig-Geo mailing list