Jonathan,

I think you can do something like this:

# remove unit from the list of arguments
# and make it a global variable instead


circ.mean.na.rm = function(x, na.rm=TRUE) {
       if(unit=="degrees")   {   x=rad(x)    }
       etc.
}

unit <-  'degrees'
azimuth_raster_agg  = aggregate(x=azimuth_raster,
fact=fact, fun=circ.mean.na.rm)

Robert

On Fri, Sep 16, 2011 at 6:21 PM, Jonathan Greenberg <jgrn@illinois.edu>wrote:

> R-sig-geo'ers:
>
> I wanted to run aggregate on a *degree* based azimuth image, with a
> "circ.mean" resampling.  Since you can't just average degrees and get
> the right answer (e.g. the arithmetic mean of 359 degrees and 1
> degrees is 180, but the "true" mean direction is 0), I adapted the
> circ.mean function in CircStats to be used with aggregate.  My
> question is, I'd like to be able to pass a "unit" argument to the
> function, but I'm not sure how to do this.  Here's the code:
>
> require(raster)
> require(CircStats)
> # require(STARStools)
> # (Swap in some DEM for the following call, or download STARStools 0.4
> or later from r-forge):
> tahoe_lidar_bareearth <-
> raster(system.file("external/tahoe_lidar_bareearth.tif",
> package="STARStools"))
> azimuth_raster <-
> slopeAspect(tahoe_lidar_bareearth,out=c('aspect'),unit='degrees')
>
> # Here is the adapted circ.mean function that adds na.rm and a "unit"
> field:
>
> circ.mean.na.rm = function(x,unit="radians",na.rm=TRUE)
> {
>        if(unit=="degrees")
>        {
>                x=rad(x)
>        }
>
>        if(na.rm)
>        {
>                x_mean=circ.mean(na.omit(x))
>        } else
>        {
>                x_mean=circ.mean(x)
>        }
>
>        if(unit=="degrees")
>        {
>                x_mean=deg(x_mean)
>        }
>        return(x_mean)
> }
>
> # How do I now pass a unit argument to aggregate?
>
> azimuth_raster_agg=aggregate(x=azimuth_raster,fact=fact,fun=circ.mean.na.rm)
> # This call returns the wrong value, because I haven't passed "unit"
> to circ.mean.na.rm() -- how do I do this from within aggregate?
>
> --j
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

