How to extract values from adjacent raster cells that are not touched by SpatialLines?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

How to extract values from adjacent raster cells that are not touched by SpatialLines?

Andre Rovai
[hidden email]

Hi all,

I've been trying to extract values from a single attribute raster (area, in
m2) that overlaps with different lines (that is, a .shp SpatialLines).

I am using the extract() function but it only extracts values from
cells that are touched (that is, crossed) by the lines.  The problem is
that my raster has adjacent cells both right and left of the cell that
is actually touched by the lines.  Thus, when I add up the extracted values
a significant amount of area (m2) is lost due to cells that were not
touched by the line and therefore values were not extracted.

I tried to work it around by:

Step 1 - first aggregating my raster to a lower resolution (i.e. increasing
the fact argument) and then
Step 2 - rasterizing the lines using this aggregated raster (created in
step 1) as a mold to make sure the rasterized lines would get thick enough
to cover the horizontal spread of cells in my original resolution raster.
Step 3 - Then I resample the rasterized lines (created in step 2) back to
the original resolution I started with.
Step 4 - Finally, extracted the values from the resampled rasterized lines
(created in step 3).

However, it didn't quite work as now the total area (m2) varies according
to the fact="" value I use when first aggregating the raster (in step 1).

I also tried using crop() and mask() but had no success as well.

I really appreciate if anyone has already dealt with a similar problem and
can help me out here.  Here are the codes I've been running to try to get
it to work:


# input raster file

g.025 <- raster("ras.asc")

g.1 <- aggregate(g.025, fact=2, fun=sum)



# input SpatialLines

Spline1 <- readOGR("/Users/xxxxx.shp")

Spline2 <- readOGR("/Users/xxxxx.shp")

Spline3 <- readOGR("/Users/xxxxx.shp")



# rasterizing using low resolution raster (aggregated)

c1 <- rasterize(Spline1, g.1, field=Spline1$type, fun=sum)

c2 <- rasterize(Spline2, g.1, field=Spline2$type, fun=sum)

c3 <- rasterize(Spline3, g.1, field=Spline3$type, fun=sum)



# resampling back to higher resolution

c1 <- resample(c1, g.025)

c2 <- resample(c2, g.025)

c3 <- resample(c3, g.025)



# preparing to extract area (m2) values from raster “g.025”

c1tab <- as.data.frame(c1, xy=T)

c2tab <- as.data.frame(c2, xy=T)

c3tab <- as.data.frame(c3, xy=T)

c1tab <- c1tab[which(is.na(c1tab$layer)!=T),]

c2tab <- c2tab[which(is.na(c2tab$layer)!=T),]

c3tab <- c3tab[which(is.na(c3tab$layer)!=T),]



# extracting area (m2) values from raster “g.025”

c1tab[,4] <- extract(g.025, c1tab[,1:2])

c2tab[,4] <- extract(g.025, c2tab[,1:2])

c3tab[,4] <- extract(g.025, c3tab[,1:2])

names(c1tab)[4] <- "area_m2"

names(c2tab)[4] <- "area_m2"

names(c3tab)[4] <- "area_m2"



# sum total area (m2)

c1_area <- sum(c1tab$area_m2)

c2_area <- sum(c2tab$area_m2)

c3_area <- sum(c3tab$area_m2)

tot_area <- sum(c1_area, c2_area, c3_area)


Thanks!

Andre Rovai
Department of Oceanography and Coastal Sciences
Louisiana State University

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: How to extract values from adjacent raster cells that are not touched by SpatialLines?

Santiago Sánchez
Hi Andre,

I'm not completely sure if this is what you are looking for, but here is a
worked example on how to get cell indices (thus, data) from neighbouring
cells. Essentially, I'm using the adjancency() function from raster,
reducing to unrepeated cells, and excluding the ones that fall under a line
object.

library(raster)
r <- raster(ncol=10,nrow=10) # generate a raster object
r[] <- 0 # populate values with 0's
p <- rasterToPolygons(r) # get a polygon for each cell (for visual purposes)

# make a line object (SpatialLines)
x <- c(-124.66110, -93.04031, -52.37858, 24.44486, 15.98469, 52.12075,
88.49592)
y <- c(-46.021148, -27.197684, -6.804443, 10.113111, 28.375197, 8.851744,
-12.463264)
xy <- data.frame(cbind(x,y))
spl <- SpatialLines(list(Lines(list(Line(xy)), ID=1)))

# get adjacent cells
lcells <- cellFromLine(r, spl)[[1]] # cells from line
r[lcells] <- 1 # mark line cells with 1
ad <- adjacent(r, lcells, 4) # this gives you a matrix with adjacent cells,
                                          # you can specify 8 or 16
neighbouring cells,
                                          # here I'm using 4
# sorting and removing duplicates
ad <- sort(as.vector(ad))
ad <- ad[!duplicated(ad)]
ad2 <- ad[ ! ad %in% lcells ]
r[ad2] <- 2

# to visualize the example:
plot(r)
plot(p, add=T)
plot(spl, add=T, col="red")

# Obviously, you con extract the values of adjacent cells simply with:
r[ad2]
[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2


Hope this helps,
Santiago

--
==========================
Santiago Sanchez-Ramirez, PhD
Postdoctoral Associate
Ecology and Evolutionary Biology
University of Toronto
==========================

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo