[R-sig-Geo] Elevation histogram for each region

[Robert J. Hijmans]

Clélia,

Here is a example with polygonValues. It works under the assumption
that all the raster cells for a single polygon can be held in memory.


library(raster)
# creating some polygons
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
pols <- SpatialPolygons( list(  Polygons(list(Polygon(p1)), 1),
Polygons(list(Polygon(p2)), 2), Polygons(list(Polygon(p3)), 3)))

# set up a raster
r <- raster(ncol=180, nrow=90)
r[] = runif(ncell(r))

# how many polygons
nzones = length(pols at polygons)

lst = list()
for (i in 1:nzones) {
	pol <- pols[i]
	p <- polygonValues(pol,r)
	lst[[i]] = hist(unlist(p), main='')  # perhaps change to
main=pol at data[1,1]  or equivalent
}

plot(lst[[1]])


Robert


On Fri, May 28, 2010 at 6:59 AM, Clélia Bilodeau
<clelia.bilodeau at gmail.com> wrote:
> Thank you for your help.
> I'm still working on it.
> I think the use of polygonValues is a good idea, but the script:
>
> library(raster)
> lst = list()
> for (i in 1:nzones) {
> p <- polygonValues(pol,r, fun=function(x){x==i})
> lst[[i]] = histogram(p[,3])
> }
> doesn't work yet.
> Do I have to define a function to use in "fun=function(x)(x==1)"?
>
> Thank you.
> C.
>
>
> 2010/5/26 Robert J. Hijmans <r.hijmans at gmail.com>:
>> Clélia,
>>
>> The optimal solution would depend a bit on your data. If the regions
>> are not too large, you could loop over the zones
>>
>> library(raster)
>> lst = list()
>> for (i in 1:nzones) {
>>  p <- rasterToPoints(r, fun=function(x){x==i})
>>  lst[[i]] = histogram(p[,3])
>> }
>>
>> If you also have the zones as polygons you could alternatively look
>> over the polygons and use  polygonValues in stead of rasterToPoints
>>
>> Another, more elaborate, option could be to, in your loop, reclassify
>> the regions raster (1 / 0 ), multiply that with the elevation data;
>> trim the output raster, make a histogram (set maxpixels to
>> ncell(raster)) (this may take a long time to finish)
>>
>> You could also first aggregate the rasters; unless you care much about
>> the extreme values.
>>
>> Hope this helps,
>> Robert
>>
>>
>> On Wed, May 26, 2010 at 2:08 AM, Clélia Bilodeau
>> <clelia.bilodeau at gmail.com> wrote:
>>> Dear list,
>>>
>>> I have two rasters, one is a Digital Elevation Model, and the other is
>>> a region file.
>>> I am looking for a method that allows me to plot the elevation
>>> histogram or density for each region.
>>>
>>> My first try was to do this with a function and to use "Zonal" from
>>> the "raster" package, but I have a very large file, so I have this
>>> error:
>>> "RasterLayers are too large. You can use fun='sum', 'mean', 'min', or
>>> 'max', but not a function"
>>>
>>> Is there another way to do this?
>>>
>>> Thank you.
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>