How to sum area of adjacent patches?

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

How to sum area of adjacent patches?

Hugo Costa
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?

Ben Tupper
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/

_______________________________________________
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?

Hugo Costa
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