[GIS] Raster-extract function – Area-Weighted values

extractr

I am working in R and want to calculate the value of a polygon. The value should consider the weights on each intersecting cell. When I try to run the "extract" function with a sample raster and polygon I get different weights that the ones I calculate manually, resulting in different final value.

Here is my sample code:

r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)    
r[] <- c(1,2,4,5)    
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)    
s.pl <- as(s, 'SpatialPolygons')    
w <- raster::extract(r, s.pl, method="simple",weights=T, normalizeWeights=F)    
mean.value <- raster::extract(r, s.pl, method="simple",weights=T, fun=mean)

The value I get is 2.14 but according to the actual weights of the cells it should be 2. Can it be because of the projection that is in lat/lon? But even when I assign projection in meters I have the same result. How can I get the value of 2 that I am interested in? I also tried with "resample" function but I get different results as well.

My final target is to create a new raster with different resolution and extents from the original one and assign the values based on the weights of the cells of the original raster that intersect with the cells of the new raster. But it seems that neither resample nor extract functions give the expected outcome.

Hope my question makes sense.

Best Answer

Based on some of the replies of this post I managed to write the following code, that works with rasterstack objects and resamples to new resolution based on the area-weighted values.

require(raster)
require(rgeos)
r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)    
r[] <- c(1,2,4,5)    
r <- stack(r, r*2, r^2)
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)    
s.pl <- as(s, 'SpatialPolygons')    
r.s <- as(r, 'SpatialPolygonsDataFrame')
pi1 <- gIntersection(r.s, s.pl, byid = T)
areas1 <- data.frame(area=sapply(pi1@polygons, FUN=function(x) {slot(x, 'area')}))
row.names(areas1) <- sapply(pi1@polygons, FUN=function(x) {slot(x, 'ID')})
areas1$Pol.old <- as.numeric(vapply(strsplit(rownames(areas1), " "), `[`, 1, FUN.VALUE=character(1)))
areas1$pol.new <- as.numeric(vapply(strsplit(rownames(areas1), " "), `[`, 2, FUN.VALUE=character(1)))
f <- r.s@data
seqs <- match(areas1$Pol.old, rownames(f))
ar <- cbind(areas1, f[seqs,])
ar[,-(1:3)] <- ar[,-(1:3)]*ar$area
f <- aggregate.data.frame(ar, by=list(ar$pol.new), FUN=sum)
f[,-(1:4)] <- f[,-(1:4)]/f$area  
ar.v <- as.matrix(f[, -c(1:4)])
s2 <- stack(s)
s1 <- setValues(s2, ar.v)

Related Solutions

[GIS] Converting data.frame with no lat-long info to raster using R

So I've created a file that reads into a data frame in the same way as yours:

> str(v2)
'data.frame':   360 obs. of  720 variables:

BUT data.frame isn't really the right thing here. Its really meant for record-oriented data, where each row is a record and each column is a potentially different variable for that record (eg each row is a person, the columns are name, age, height, etc).

So you really only need to scan the data in as one long vector and feed it to a raster.

Step 1, define an empty raster of the right size and shape (note I'm assuming the raster covers the whole world, so the limits are not the cell centres):

> m2=raster(nrow=360,ncol=720,xmn=-180,xmx=180,ymn=-90,ymx=90)

Step 2, read numeric values into the raster data slot:

> m2[]=scan("d.txt",what=1)
Read 259200 items

And give it a projection if needed:

> projection(m2)="+init=epsg:4326"
> plot(m2)

If you want to check that the resolution and the cell centres are as expected, use these functions:

> res(m2)
[1] 0.5 0.5
> xFromCol(m2,1:10)
 [1] -179.75 -179.25 -178.75 -178.25 -177.75 -177.25 -176.75 -176.25 -175.75
[10] -175.25
> yFromRow(m2,1:10)
 [1] 89.75 89.25 88.75 88.25 87.75 87.25 86.75 86.25 85.75 85.25

which shows the resolution is half a degree and the cell centres (or at least the first 10) are at those specified coordinates.

[GIS] Extract the edge between two raster cells with different values

I think it is not wise to mix the clump(..) functionality from igraph with the dissolve=TRUE parameter from the rasterToPolygon routine. They both do something with to aggregate the fields together but in a different way. At least we want to do 3 things:

read or desing a raster
select raster area where the contour goes around
define how the contour is shaped (rectangular,linear or smooth) and define which areas belong together.

In the clump code raster, the type of raster and the definition of NA seems to be important to steer clumping the process. I made some test but with bad results. I followed your sketch and here is a little analysis how to get things work:

# Load packages
require('raster')
require('rgeos')

# Clean up everything
rm(list=ls())

# Set a defined random seed
set.seed(2)

# Create a float raster with values
# the interval [0,2] float
rs <- raster(nrow=10, ncol=10)

# Scale the random numbers from interval [0,1) to
# to [0,2.2) and shift the interval to [-0.1,2.1)
rs[] <- runif(ncell(rs)) * 2.2 - 0.1

# Cut off the raster values to [0,2]
# Everything smaller then 0 is zero
values(rs)[values(rs) < 0.0] <- 0.0

# Everything larger then 2 is two
values(rs)[values(rs) > 2.0] <- 2.0

Is the raster field well constructed?

# Construction is OK?
> quantile(values(rs))
       0%       25%       50%       75%      100% 
0.0000000 0.3928834 0.8733250 1.6027794 2.0000000

Here is the function prototype that selects the field in the interval [1,2).

# Function of the contour ID
inOne <- function(x) { x>=1 & x<2 }

> inOne(0) 
[1] FALSE

> inOne(1) 
[1] TRUE

> inOne(2)
[1] FALSE

The contour cannot be dissolved (clump) because of the float number nature of the raster field is distinct in the ID process (dissolve-=TRUE).

# Contour of the float desing
# x := [1,2)
ct <- rasterToPolygons(rs, 
         fun=inOne , 
         dissolve=TRUE) 
plot(rs)
plot(ct, add=TRUE)

You see the right contour groups, but polygons are not joint.

So if we have a unique ID of each cell as Integer, the dissolve process should work.

# Apply integer operation (ceiling, floor, round)
# to the float number fields 
rs.int <- ceiling(rs)
values(rs.int)
ct.int <- rasterToPolygons(rs.int, 
              fun=inOne , 
              dissolve=TRUE) 
plot(rs.int)
plot(ct.int, add=TRUE)

Conclusion: I think (do know not exactly) what behind the clump stuff works a raster based region growing routine. The process dissolve=TRUE in the rasterToPolygon (based on rgeos CRAN) seems to follow a vector approach. So I've to read the manuals of igraph and rgeos carfully.

REM: The selection of contours (float vs. int) differs, because of the nature of (ceiling, floor and round).

Best Answer

Related Solutions

[GIS] Converting data.frame with no lat-long info to raster using R

[GIS] Extract the edge between two raster cells with different values

Related Question