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

[Robert J. Hijmans]

Can you avoid transforming, by computing the surface area for each
lon/lat cell, and then extract the values by polygon. Perhaps
something like this (I am not sure if this is how slope correction
should be done)?

library(raster)
elevation <- getData('alt', country='CHE')
slope <- terrain(elevation, opt='slope', unit='radians')
# c^2 = a^2 + b^2; let a be 1; atan(slope) = b/a = b
correction <- sqrt(atan(slope)^2 + 1)
area <- area(elevation)
cor_area <- area * correction

Robert

On Sun, Dec 22, 2013 at 1:38 PM, Michael Sumner <mdsumner at gmail.com> wrote:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo