Moving-window sum in a matrix

4 messages
Open this post in threaded view
|

Moving-window sum in a matrix

 Dear colleagues, I have a 500 by 500 matrix with presence or absence value in each cell as an input matrix. I am trying to create an output matrix having in each cell the sum of the values from neigbouring cells in the input matrix according to a moving-window defined as: 0 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 0 where the output cell receiving the values is located in the center. I have written a function that loop over the cells of the output matrix, and sum up the values from these neigbouring cells. It works, but it is fairly slow, and hangs if I try it over a larger matrix. Given that this is a rather basic function available in most raster GIS packages (IDRISI and AV Spatial Analyst) and running fairly quickly, I wondered whether anyone had previously come accross a similar function in an existing R-package. Any helps  or pointer would be appreciated, Many thanks in advance, Marius
Open this post in threaded view
|

Moving-window sum in a matrix

 Hello, I think that you could do it using spdep package. I suppose that, if you: - construct a neighborhoord object (nb class). You can define a distance based neighborhood (using dnearneigh) as it seems that your definition is similar to say i and j neighbors if 0 myvalues=rbinom(100,size=1,prob=0.5)  > xy=expand.grid(1:10,1:10)  > dnn1=dnearneigh(as.matrix(xy),0,sqrt(1+4))  > plot(dnn1,as.matrix(xy))  > nb2listw(dnn1,style="B")  > lw1=nb2listw(dnn1,style="B")  > sum.neigh=lag.listw(lw1,myvalues)  > sum.neigh    [1]  3  6  4  5  7  5  4  4  3  3  6  8  6  7  9  8  6  8  6  4  5  9 10 11 10  7  8  9  7  5   [31]  7 10 14 11 10  8  8  9  8  5  8 12 11 10 11 12 10 12 11  7  9 11 12  9 13 11 10 12 12  9   [61]  6 11 12 13 10 12 10 13  9  7  6  9 12 10 12 10 13 12  9  5  6  8  8 10 10 11  8  9  7  5   [91]  4  6  7  6  9  7  6  4  5  3 Roger could confirm that it is correct ! At 10:05 16/12/2004, Marius Gilbert wrote: >Dear colleagues, > >I have a 500 by 500 matrix with presence or absence value in each cell as >an input matrix. I am trying to create an output matrix having in each >cell the sum of the values from neigbouring cells in the input matrix >according to a moving-window defined as: > >0 1 1 1 0 >1 1 1 1 1 >1 1 0 1 1 >1 1 1 1 1 >0 1 1 1 0 > >where the output cell receiving the values is located in the center. I >have written a function that loop over the cells of the output matrix, and >sum up the values from these neigbouring cells. It works, but it is fairly >slow, and hangs if I try it over a larger matrix. Given that this is a >rather basic function available in most raster GIS packages (IDRISI and AV >Spatial Analyst) and running fairly quickly, I wondered whether anyone had >previously come accross a similar function in an existing R-package. > >Any helps  or pointer would be appreciated, > >Many thanks in advance, > >Marius > >_______________________________________________ >R-sig-Geo mailing list >R-sig-Geo at stat.math.ethz.ch >https://stat.ethz.ch/mailman/listinfo/r-sig-geoSt?phane DRAY -------------------------------------------------------------------------------------------------- D?partement des Sciences Biologiques Universit? de Montr?al, C.P. 6128, succursale centre-ville Montr?al, Qu?bec H3C 3J7, Canada Tel : (514) 343-6111 poste 1233         Fax : (514) 343-2293 E-mail : stephane.dray at umontreal.ca -------------------------------------------------------------------------------------------------- Web                                          http://www.steph280.freesurf.fr/
Open this post in threaded view
|

Moving-window sum in a matrix

 Stephane, have you tried how well this scales to, say, 1000x1000 matrices? Indeed, most basic raster GIS functionality is not available in R packages, and would be welcome. R code is efficient for small problems, but for cases like this scales badly as Marius noted. We are actually working on an interface from R to the functionality of PCRaster (http://pcraster.geog.uu.nl -- yet another raster GIS), and expect to have this working within a few months. PCRaster is not open source, so it is unlikely to become a CRAN package straight away, but still it may fill a need. This all depends on availability of sp, the long awaited package with S classes for spatial point, line, grid, arc/node and polygon data. We're working very hard on this, and most likely will announce a beta in January. I expect a PCRasteR package to follow that soon. -- Edzer Stephane DRAY wrote: > Hello, > I think that you could do it using spdep package. I suppose that, if you: > - construct a neighborhoord object (nb class). You can define a distance > based neighborhood (using dnearneigh) as it seems that your definition > is similar to say i and j neighbors if > 0 - then transform the nb to listw with nb2listw with style="B" > - then use lag.listw to sompute the sum > > an example: >  > myvalues=rbinom(100,size=1,prob=0.5) >  > xy=expand.grid(1:10,1:10) >  > dnn1=dnearneigh(as.matrix(xy),0,sqrt(1+4)) >  > plot(dnn1,as.matrix(xy)) >  > nb2listw(dnn1,style="B") >  > lw1=nb2listw(dnn1,style="B") >  > sum.neigh=lag.listw(lw1,myvalues) >  > sum.neigh >   [1]  3  6  4  5  7  5  4  4  3  3  6  8  6  7  9  8  6  8  6  4  5  9 > 10 11 10  7  8  9  7  5 >  [31]  7 10 14 11 10  8  8  9  8  5  8 12 11 10 11 12 10 12 11  7  9 11 > 12  9 13 11 10 12 12  9 >  [61]  6 11 12 13 10 12 10 13  9  7  6  9 12 10 12 10 13 12  9  5  6  8   > 8 10 10 11  8  9  7  5 >  [91]  4  6  7  6  9  7  6  4  5  3 > > Roger could confirm that it is correct ! > > > > At 10:05 16/12/2004, Marius Gilbert wrote: > >> Dear colleagues, >> >> I have a 500 by 500 matrix with presence or absence value in each cell as >> an input matrix. I am trying to create an output matrix having in each >> cell the sum of the values from neigbouring cells in the input matrix >> according to a moving-window defined as: >> >> 0 1 1 1 0 >> 1 1 1 1 1 >> 1 1 0 1 1 >> 1 1 1 1 1 >> 0 1 1 1 0 >> >> where the output cell receiving the values is located in the center. I >> have written a function that loop over the cells of the output matrix, >> and >> sum up the values from these neigbouring cells. It works, but it is >> fairly >> slow, and hangs if I try it over a larger matrix. Given that this is a >> rather basic function available in most raster GIS packages (IDRISI >> and AV >> Spatial Analyst) and running fairly quickly, I wondered whether anyone >> had >> previously come accross a similar function in an existing R-package. >> >> Any helps  or pointer would be appreciated, >> >> Many thanks in advance, >> >> Marius >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo at stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo> > > St?phane DRAY > -------------------------------------------------------------------------------------------------- > > D?partement des Sciences Biologiques > Universit? de Montr?al, C.P. 6128, succursale centre-ville > Montr?al, Qu?bec H3C 3J7, Canada > > Tel : (514) 343-6111 poste 1233         Fax : (514) 343-2293 > E-mail : stephane.dray at umontreal.ca > -------------------------------------------------------------------------------------------------- > > Web                                           > http://www.steph280.freesurf.fr/> > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo at stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Open this post in threaded view
|

Hex GIS, (was) Moving-window ...

 Does anyone know of a GIS based on a regular hexagon tessellation rather than on squares or rectangles?  (There was a project twenty some years ago called System X developed for the Navy by Dean Lucas that started into commercialization but never finished as far as I know.) r-sig-geo-bounces at stat.math.ethz.ch wrote on 2004-12-16 12:17:28: > Stephane, have you tried how well this scales to, say, 1000x1000 > matrices? Indeed, most basic raster GIS functionality is not available > in R packages, and would be welcome. R code is efficient for small > problems, but for cases like this scales badly as Marius noted. > > We are actually working on an interface from R to the functionality of > PCRaster (http://pcraster.geog.uu.nl -- yet another raster GIS), and > expect to have this working within a few months. PCRaster is not open > source, so it is unlikely to become a CRAN package straight away, but > still it may fill a need. > > This all depends on availability of sp, the long awaited package with S > classes for spatial point, line, grid, arc/node and polygon data. We're > working very hard on this, and most likely will announce a beta in > January. I expect a PCRasteR package to follow that soon. > -- > Edzer