[R-sig-Geo] Raster - zonal statistics -zones are not matching the mask?

[Robert J. Hijmans]

Byman,

I suspect that raster cell sizes are large relative to some of the
polygons (small island states), and that not all polygons cover at
least one raster cell. You can check that with:

um <- unique(msk)
length(um)

# what you are doing can be done more directly, I think, like this:

library(raster)
world <- readShapePoly("world.shp", proj4string=CRS("+proj=longlat"))
bdat <- brick("global.obs.Prcp.1950-1999.nc",varname="Prcp")

msk <- rasterize(world, bdat)
zst <- zonal(bdat, msk, stat='median')

# and now perhaps this?
zst[,2:ncol(zst)] <- zst[,2:ncol(zst)] * 12 * 30


# to deal with the missing values, you can do something like this
x <- data.frame(zone=1:480)
zst2 <- merge(x, zst, by='zone', all.x=T)

world at data <- cbind(world at data, zst2[,-1])

# To insist in getting a value for the small countries, you can do
this, with a recent version of raster
# (but this is much slower than zonal)

x <- extract(bdat, world, fun=median, small=TRUE)

Best, Robert


On Sun, Jun 5, 2011 at 3:22 AM, Byman <bymanh at gmail.com> wrote:
> Hi
>
> I am computing average values on country /state basis globally. I am using
> raster package in R to compute the zonal values &  I read a shape file from
> ArcGIS in R also. Here is my short process i used. I am not able to get the
> zones to match the polygon in the mask.....
>
>     world<-readShapePoly("world.shp", proj4string=CRS("+proj=longlat"))
>
> bdat<-brick("global.obs.Prcp.1950-1999.nc",varname="Prcp")
> # read the NetCDF file uisng raster package
>     ras<-raster(ncol=bdat at ncols,nrow=bdat at nrows)
>                     #create mask with same attributes with raster
>     ras[]<-1
>     extent(ras)<-extent(world)
>     msk <- rasterize(world, ras)
>     bdat <- mask(bdat, msk)
>
> zst<-zonal(bdat*12*30,msk,stat='median')
> #Compute the zonal statistics
>
> The shape file has 480 polygons and I expect to get the same number of zones
> but to my surprise I get 337 zones..  As such I cannot map back/assign the
> results of the zonal function to the object in order to plot. The mask (msk)
> has the 480 zones but after the step I get different number of zones. What
> am I missing here.
>
> world$zonal = zst
>
> Error in `[[<-.data.frame`(`*tmp*`, name, value = list(zone = c(3, 4,  :
>   replacement has 337 rows, data has 480
>
> I have not supplied the files but I can do so.
>
> Thank you in advance for suggested solutions
>
> rgds
> Byman
>