[R-sig-Geo] Calculate (global/local) Moran I for a segmented satellite image

[Ruben Van De Kerchove]

Dear,

I am trying to calculate global Moran I and local Moran I for a
segmented satellite image.  This is what I would like to obtain.

- Segment a high resolution satellite in GRASS using i.segment.
- Calculate mean band values for each created object
- Calculate local (for each zone) and global Moran I for these mean
object values taking into account only the neighboring objects

Step 1 works fine, but I have some issues with steps 2 and 3.
Therefore, below you find an example I created for these steps.

Although the code seems to work on this example, it isn’t optimized in
terms of speed (and memory?) for very large raster sets. Especially
rasterToPolygons() is extremely slow. I think that it should be better
to calculate neighbors immediately on the “raster regions”, but is
this possible? I checked cell2nb() but this seems not to work for
regions/objects. Furthermore it seems that when I convert large
rasters (100 000 000) cells into polygons the number of created
polygons might slightly differ from the number of unique regions in
the raster (tested with gdal_polygonize and r.to.vect; not yet with
rasterToPolygons due to the time it takes). Any idea why this happens?

Would it be possible to use Moran and MoranLocal from the raster
package in this case (maybe in combination with the zonal function? I
guess that this would also speed up things.

Thanks and regards,

Ruben

----------------------------------------------
library(spgrass6)
library(raster)
library(spdep)

###LOAD DATA

#Satellite image
r <- c(55, 56, 54, 58, 65, 75, 82, 77, 74, 74, 69, 61, 62, 71, 73,
         63, 62, 63, 64, 59, 85, 88, 95, 106, 110, 99, 89, 82, 79, 84, 97, 79,
         55, 55, 56, 60, 91, 95, 86, 98, 115, 105, 110, 107, 101, 89, 85, 68,
         55, 54, 53, 61, 82, 102, 88, 93, 96, 94, 110, 114, 109, 103, 92, 68,
         59, 58, 60, 64, 88, 99, 82, 81, 71, 80, 89, 89, 89, 102, 104, 75, 63,
         57, 58, 58, 77, 92, 82, 71, 59, 90, 105, 92, 79, 98, 110, 83, 62, 55,
         56, 56, 80, 90, 99, 88, 64, 91, 112, 94, 76, 91, 100, 85, 62, 59, 55,
         61, 99, 97, 93, 80, 65, 87, 107, 80, 59, 60, 67, 66, 65, 62, 68, 72,
         102, 94, 90, 83, 74, 81, 96, 69, 52, 50, 51, 54, 62, 62, 86, 85, 55,
         59, 64, 72, 75, 70, 70, 62, 66, 61, 55, 57, 52, 59, 61, 56, 41, 40,
         43, 44, 46, 48, 50, 52, 68, 69, 60, 61, 42, 43, 44, 43, 42, 41, 42,
         44, 43, 43, 44, 47, 58, 59, 55, 61, 44, 41, 39, 42, 44, 43, 42, 42,
         42, 43, 43, 49, 56, 53, 53, 61, 43, 42, 40, 42, 42, 42, 41, 42, 42,
         43, 42, 53, 66, 61, 51, 62, 40, 42, 41, 40, 43, 49, 46, 42, 42, 43,
         43, 49, 59, 62, 53, 62, 40, 41, 42, 43, 49, 54, 47, 44, 42, 44, 43,
         46, 52, 56, 56, 61)
r <- matrix(r, ncol=16, nrow=16, byrow=T)
r<- raster(r)
plot(r)

#Segments
seg <- r
k=0
for (i in 1:4){
  for (j in 1:4){
    k=k+1
  seg[(4*i-3):(4*i),(4*j-3):(4*j)]=matrix(rep(k,16),4,4)
}
}
plot(seg)

####ZONAL STAT
#Create raster of mean for each segment
r.mean=zonal(r,seg, fun="mean")
r.mean=as.data.frame(r.mean)
r.mean.im=subs(seg,r.mean,by="zone",which="mean")
plot(r.mean.im)

#Create polygons of mean for each segment
r.mean.pl=rasterToPolygons(seg,dissolve=T)
r.mean.pl=spCbind(r.mean.pl,r.mean[,2])
names(r.mean.pl)[2]="mean"

#####MI
#Calculate neighbours and weights
nb=poly2nb(r.mean.pl)
plot(r.mean.pl, border="grey")
plot(nb, coordinates(r.mean.pl), add=TRUE)
nb.weight=nb2listw(nb,style="B")

values=as.data.frame(r.mean.pl)[,2]
MI.global=moran(values,nb.weight,length(nb),Szero(nb.weight))
MI.local=localmoran(values,nb.weight)

r.mean.pl=spCbind(r.mean.pl,as.data.frame(MI.local)[[1]])
names(r.mean.pl)[3]="MI.local"

spplot(r.mean.pl, zcol="MI.local", col=NA)