# How to sum area of adjacent patches? Classic List Threaded 3 messages Reply | Threaded
Open this post in threaded view
|

## How to sum area of adjacent patches?

 Deal all, I'm merging several rasters representing patches of forest, and the pixel values represent the area of the patches. The patches' area should be summed along the edges of the merged rasters. I provide a reproducible example. Any help is welcome. Hugo library(raster) r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1) r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2) r1[]<-c(rep(24,24), rep(0,76)) r2[]<-c(rep(0,86), rep(14,14)) par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2") r3<-merge(r1,r2) ; plot(r3, main="merged rasters") # do something r3[r3>0]<-38     ; plot(r3, main="desired result")         [[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 sum area of adjacent patches?

 Hi, I think you may need to use raster::clump() and table() to label your regions.  Check out an old message about clump... https://stat.ethz.ch/pipermail/r-sig-geo/2017-January/025346.htmlIt's sort of cumbersome, but I think the following works... library(raster) r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1) r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2) r1[]<-c(rep(24,24), rep(0,76)) r2[]<-c(rep(0,86), rep(14,14)) par(mfrow=c(3,2)); plot(r1,main="r1"); plot(r2,main="r2") r3<-merge(r1,r2) ; plot(r3, main="merged rasters") r3 <- r3 r4[r4>0]<-38     ; plot(r4, main="desired result") r5 <- clump(r3) ; plot(r5, main = 'clumped') r5[is.na(r5)] <- 0 r5t <- table(getValues(r5)) for (i in names(r5t)[-1]){         id <- as.numeric(i)         print(id)         r5[r5 == id] <- r5t[i] } plot(r5, main = 'relabeled with area') Cheers, Ben P.S. Years ago I used a language called IDL.  It's histogram function has an output keyword called reverse_indices (more here https://www.harrisgeospatial.com/docs/HISTOGRAM.html#H_835179117_677188 and here http://www.idlcoyote.com/tips/histogram_tutorial.html )  - super useful!  Your situation just begs for that same solution in R. > On Nov 23, 2018, at 5:13 PM, Hugo Costa <[hidden email]> wrote: > > Deal all, > > I'm merging several rasters representing patches of forest, and the pixel > values represent the area of the patches. The patches' area should be > summed along the edges of the merged rasters. I provide a reproducible > example. Any help is welcome. > Hugo > > library(raster) > r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1) > r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2) > r1[]<-c(rep(24,24), rep(0,76)) > r2[]<-c(rep(0,86), rep(14,14)) > par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2") > > > r3<-merge(r1,r2) ; plot(r3, main="merged rasters") > > # do something > > r3[r3>0]<-38     ; plot(r3, main="desired result") > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geoBen Tupper Bigelow Laboratory for Ocean Sciences 60 Bigelow Drive, P.O. Box 380 East Boothbay, Maine 04544 http://www.bigelow.orgEcological Forecasting: https://eco.bigelow.org/_______________________________________________ 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 sum area of adjacent patches?

 Hi Ben, thank you for your help. Clump is a possibility, however I was looking for something different. In my real data, 9 rasters have to be merged (1 central raster and the 8 surrounding rasters), and running clump after merging 9 rasters is slow and memory consuming. And this has to be repeated thousands of times (hopefully in parallel with the rasters in memory). I didn't explain these details to make the message short. I'll continue searching for more alternatives, and all ideas are welcome. Thanks Hugo Ben Tupper <[hidden email]> escreveu no dia sexta, 23/11/2018 à(s) 23:58: > Hi, > > I think you may need to use raster::clump() and table() to label your > regions.  Check out an old message about clump... > > https://stat.ethz.ch/pipermail/r-sig-geo/2017-January/025346.html> > It's sort of cumbersome, but I think the following works... > > > library(raster) > r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1) > r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2) > r1[]<-c(rep(24,24), rep(0,76)) > r2[]<-c(rep(0,86), rep(14,14)) > par(mfrow=c(3,2)); plot(r1,main="r1"); plot(r2,main="r2") > > > r3<-merge(r1,r2) ; plot(r3, main="merged rasters") > > r3 <- r3 > r4[r4>0]<-38     ; plot(r4, main="desired result") > > r5 <- clump(r3) ; plot(r5, main = 'clumped') > r5[is.na(r5)] <- 0 > > r5t <- table(getValues(r5)) > for (i in names(r5t)[-1]){ >         id <- as.numeric(i) >         print(id) >         r5[r5 == id] <- r5t[i] > } > > plot(r5, main = 'relabeled with area') > > > Cheers, > Ben > > P.S. Years ago I used a language called IDL.  It's histogram function has > an output keyword called reverse_indices (more here > https://www.harrisgeospatial.com/docs/HISTOGRAM.html#H_835179117_677188> and here http://www.idlcoyote.com/tips/histogram_tutorial.html )  - super > useful!  Your situation just begs for that same solution in R. > > > > > > > On Nov 23, 2018, at 5:13 PM, Hugo Costa <[hidden email]> wrote: > > > > Deal all, > > > > I'm merging several rasters representing patches of forest, and the pixel > > values represent the area of the patches. The patches' area should be > > summed along the edges of the merged rasters. I provide a reproducible > > example. Any help is welcome. > > Hugo > > > > library(raster) > > r1<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=0, ymx=1) > > r2<-raster(nrows=10, ncols=10, xmn=0, xmx=1, ymn=1, ymx=2) > > r1[]<-c(rep(24,24), rep(0,76)) > > r2[]<-c(rep(0,86), rep(14,14)) > > par(mfrow=c(2,2)); plot(r1,main="r1"); plot(r2,main="r2") > > > > > > r3<-merge(r1,r2) ; plot(r3, main="merged rasters") > > > > # do something > > > > r3[r3>0]<-38     ; plot(r3, main="desired result") > > > >       [[alternative HTML version deleted]] > > > > _______________________________________________ > > R-sig-Geo mailing list > > [hidden email] > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo> > Ben Tupper > Bigelow Laboratory for Ocean Sciences > 60 Bigelow Drive, P.O. Box 380 > East Boothbay, Maine 04544 > http://www.bigelow.org> > Ecological Forecasting: https://eco.bigelow.org/> > > > > >         [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo