[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 |
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 |
Free forum by Nabble | Edit this page |