[R-sig-Geo] Mosaic Option to BLEND

[Rafael Wüest]

Dear Etienne and Robert

building on Etienne's idea, I was just about to get to Robert's suggestion.

Thank you very much, I will do some further testing, but this is pretty much what I was looking for. Thanks to both of you.

Best,
Rafael

On 25.10.2012, at 19:19, Robert J. Hijmans wrote:

> 
> Some ideas, building on Etienne's script:
> 
> 
> library(raster)
> m <- matrix(NA, 50, 50)
> r1 <- raster(m)
> r1[] <- 10
> r2 <- shift(r1, x = 0.4, y = 0.4)
> r2[] <- 20
> 
> 
> i <- intersect(raster(r1), raster(r2))
> j <- raster::expand(i, c(1,1)) 
> a <- crop(r1, j)
> values(a) <- 1
> b <- crop(r2, j)
> values(b) <- 2
> ab <- merge(a, b)
> ba <- merge(b, a)
> p1 <- rasterToPoints(ab, function(x) x==2)
> p2 <- rasterToPoints(ba, function(x) x==1)
> 
> xy <- xyFromCell(i, 1:ncell(i))
> d1 <- distanceFromPoints(i, p1[,1:2])
> d2 <- distanceFromPoints(i, p2[,1:2])
> dsum <- d1 + d2
> 
> z1 <- d1 * crop(r1, d1) / dsum
> z2 <- d2 * crop(r2, d2) / dsum
> m <- merge(z1 + z2, r1, r2)
> 
> plot(m)
> 
> 
> 
> 
> 
> On Thu, Oct 25, 2012 at 8:40 AM, Etienne B. Racine <etiennebr at gmail.com> wrote:
> I find it not very obvious what ESRI BLEND option is doing.
> 
>    - BLEND —The output cell value of the overlapping areas will be a
>    horizontally weighted calculation of the values of the cells in the
>    overlapping area. [
>    http://help.arcgis.com/fr/arcgisdesktop/10.0/help/index.html#//00pv00000007000000
>    ]
> 
> or 9.2
> 
> 
>    - BLEND — The output cell value of the overlapping areas will be a blend
>    of values that overlap; this blend value is based on an algorithm that is
>    weight based and is dependent on the distance from the pixel to the edge
>    within the overlapping area. [
>    http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=mosaic_to_new_raster_%28data_management%29
>    ]
> 
> I think everything is about how you generate the weight matrix here.  I
> couldn't find with a quick search, what horizontal weight is. But let's say
> it is about giving a "priority" to cells closer to their center this makes
> sense if the overlap of the rasters is less than 50 % in each axis.
> # generate data
> library(raster)
> m <- matrix(NA, 5, 5)
> r <- raster(m)
> r[] <- 1
> r2 <- shift(raster(m), x = 0.4, y = 0.4)
> r2[] <- 2
> 
> # create a distance-weighted matrix
> d <- raster(m)
> d2 <- r2
> 
> # This is a simple example
> # Everything is happening here. How would you generate a weight
> # that make sense for your context ?
> d[] <- 1 / (as.matrix(dist(xyFromCell(r, 1:ncell(r))))[13, ] + 1e-9)
> d2[] <- 1 / (as.matrix(dist(xyFromCell(r2, 1:ncell(r))))[13, ] + 1e-9)
> 
> # merge weights
> w <- mosaic(d, d2, fun = sum)
> # create weighted rasters
> wr <- r * d / w
> wr2 <- r2 * d2 / w
> # simply sum the two weighted rasters
> mos <- mosaic(wr, wr2, fun = sum)
> plot(mos)
> 
> However as I mentioned earlier it is not working very well if the rasters
> overlap more than 50% in an axis or the other since the distance from the
> center start to decrease. Try it with a smaller shift, e.g. r2 <-
> shift(raster(m), x = 0.2, y = 0.2).
> 
> Maybe it will get you started however
> 
> Etienne
> 
> 2012/10/25 Rafael Wüest <rafael.wueest at wsl.ch>
> 
> > Dear R geo experts,
> >
> > I have specific question to the raster::mosaic function.
> >
> > From the help file of the mosaic function (raster version 2.0-08), I see
> > that we have options to assign the min, max, mean to the mosaic'ed cells.
> > From ArcINFO, I'm used to the BLEND option, which will assign a value to
> > the mosaic'ed cell that relies on an algorithm that is weight based and
> > dependent on the distance from the pixel to the edge within the overlapping
> > area.
> >
> > Now, is there an easy way I can get the same result using the
> > raster::mosaic function? Or is there another elegant way to achieve the
> > same using some other packages or functions?
> >
> > I'm grateful for any hints!
> >
> > Rafael
> >
> >
> > --
> > Rafael Wüest
> > Swiss Federal Research Institute WSL
> > Zürcherstrasse 111
> > 8903 Birmensdorf
> > Switzerland
> >
> > +41 44 7392126
> > rafael.wueest at wsl.ch
> > http://www.wsl.ch/info/mitarbeitende/wueest/index_EN
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> 
>         [[alternative HTML version deleted]]
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
> 



--
Rafael Wüest
Swiss Federal Research Institute WSL
Zürcherstrasse 111
8903 Birmensdorf
Switzerland

+41 44 7392126
rafael.wueest at wsl.ch
http://www.wsl.ch/info/mitarbeitende/wueest/index_EN