[R-sig-Geo] Map of pixel area?

[Jonathan Greenberg]

Thanks all!  I ended up doing more along the lines of what Forrest
recommended (the rectangular assumption, for a global analysis, seemed
inappropriate).  Here's my solution using rasterEngine (spatial.tools) that
converts to square kilometers:

pixel_area_calculator <- function(x,CRS="+proj=moll
+ellps=WGS84",scaling=1/1000000)
{
x <- raster(x)
 pixel_area_function <- function(x,CRS,scaling)
{
x_poly <- rasterToPolygons(x,n=16)
x_tr <- spTransform(x_poly,CRS=CRS(CRS))
x_tr_area <- gArea(x_tr,byid=TRUE)
dim(x_tr_area) <- c(ncol(x),nrow(x),1)
x_tr_area <- x_tr_area*scaling
return(x_tr_area)
}
 pixel_area <-
rasterEngine(x=x,fun=pixel_area_function,debugmode=FALSE,chunksize=1,
chunk_format="raster",args=list(CRS=CRS,scaling=scaling),.packages="rgeos")
return(pixel_area)
}

On Fri Jan 09 2015 at 5:21:34 PM Michael Sumner <mdsumner at gmail.com> wrote:

> See ?area in raster for the standard approximation method.
>
> Cheers, Mike
>
> On Sat, 10 Jan 2015 10:15 Forrest Stevens <forrest at ufl.edu> wrote:
>
>> Oops, sorry for the mistake when copying my code in.  Please find the
>> working code for the areas here:
>>
>> r1 <- raster()
>> r1[] <- 1:ncell(r1)
>> p1 <- rasterToPolygons(r1)
>>
>> ## Project to cylindrical equal area:
>> require(rgdal)
>> p2 <- spTransform(p1, CRS("+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30
>> +a=6371228.0 +units=m"))
>>
>> ## Extract areas and assign back to raster:
>> areas <- sapply(slot(p2, "polygons"), slot, "area")
>> r2 <- r1
>> r2[] <- areas
>> plot(r2)
>>
>>
>> Sorry for the confusion!
>>
>> Sincerely,
>> Forrest
>>
>> --
>> Forrest R. Stevens
>> Ph.D. Candidate, QSE3 IGERT Fellow
>> Department of Geography
>> Land Use and Environmental Change Institute
>> University of Florida
>> www.clas.ufl.edu/users/forrest
>>
>>
>> On Fri, Jan 9, 2015 at 6:11 PM, Forrest Stevens <forrest at ufl.edu> wrote:
>> > The easiest way I can think of would be to create an "index" raster
>> > based on your raster of interest, convert it to polygons, project it
>> > to an equal area projection with linear units amenable to you, then
>> > extract the polygon areas and assign them back to your index raster
>> > for weighting.  Perhaps something like this:
>> >
>> >
>> > r1 <- raster()
>> > r1[] <- 1:ncell(r1)
>> > p1 <- rasterToPolygons(r1)
>> >
>> > ## Project to cylindrical equal area:
>> > require(rgdal)
>> > p2 <- spTransform(p1, CRS("+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30
>> > +a=6371228.0 +units=m"))
>> >
>> > ## Extract areas and assign back to raster:
>> > areas <- sapply(slot(p2, "polygons"), slot, "area")
>> > r2 <- r1
>> > plot(r2)
>> >
>> >
>> > Though, maybe someone else has a better way?
>> >
>> > Sincerely,
>> > Forrest
>> > --
>> > Forrest R. Stevens
>> > Ph.D. Candidate, QSE3 IGERT Fellow
>> > Department of Geography
>> > Land Use and Environmental Change Institute
>> > University of Florida
>> > www.clas.ufl.edu/users/forrest
>> >
>> >
>> > On Fri, Jan 9, 2015 at 5:29 PM, Jonathan Greenberg <jgrn at illinois.edu>
>> wrote:
>> >> R-sig-geo'ers:
>> >>
>> >> I have an analysis where I need to take a geographic projection map,
>> e.g.:
>> >>
>> >> r1 <- raster()
>> >>
>> >> And calculate, for each pixel, the square kilometers that each
>> individual
>> >> pixel covers (since they vary across a geographic projection) -- the
>> idea
>> >> is to end up with a map of pixel area.  Can anyone think of a clever
>> way to
>> >> do this?  I'll be using this to do some area weighted summaries, but I
>> >> can't simply reproject the input map.  Thanks!
>> >>
>> >> --j
>> >>
>> >>         [[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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>

	[[alternative HTML version deleted]]