Moving-window sum in a matrix

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Moving-window sum in a matrix

Marius Gilbert
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



Reply | Threaded
Open this post in threaded view
|

Moving-window sum in a matrix

Stéphane Dray
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<dij<sqrt(cellsize^2+(2*cellsize)^2)
- 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/



Reply | Threaded
Open this post in threaded view
|

Moving-window sum in a matrix

Edzer J. Pebesma
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<dij<sqrt(cellsize^2+(2*cellsize)^2)
> - 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



Reply | Threaded
Open this post in threaded view
|

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

White.Denis@epamail.epa.gov
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