[R-sig-Geo] Memory problems with dnearneigh in spdep

Wed Apr 10 09:10:41 CEST 2013

On Tue, 9 Apr 2013, Empty Empty wrote:

> Thanks so much, Roger, for your suggestions.
> I should have realized it was taking kilometers rather than meters!
>
> My summary(r) is:
> Object of class SpatialPointsDataFrame
> Coordinates:
>                 min      max
> coords.x1 -17.52090 51.40407
> coords.x2 -34.82918 37.53746
> Is projected: TRUE 
> proj4string :
> [+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84
> +datum=WGS84 +units=m +no_defs +towgs84=0,0,0]

Given your object summary, you need to check where your data came from. 
The bounding box is for geographical coordinates, but the declared 
coordinate reference system is projected in units of metres. So your 9.333 
is 9.333m, and no neighbours will be found. The inclusion of NAs in 
gridded data represented as points is unnecessary, and all points not on 
land should be dropped before analysis begins.

I do suggest moving to a raster representation, and using focal methods in 
the raster package to generate the separate components of equation 14.3, 
p. 263-264 in Gettis & Ord (1996). Using focal methods defines the moving 
window as a matrix, here 10x10, moved over the raster and circumventing 
the creation of a weights object - create a new raster layer with \sum{j} 
w_{ij}(d) x_j values. W_i^* will be a constant, as will \bar{x}^*, and 
possibly the other starred terms too.

Hope this helps,

Roger

> Number of points: 36820887
> Data attributes:
>     POINTID           GRID_CODE        
>  Min.   :      1    Min.   :     0.00  
>  1st Qu.: 250000    1st Qu.:     0.28  
>  Median : 500000    Median :     2.57  
>  Mean   : 500000    Mean   :    26.43  
>  3rd Qu.: 750000    3rd Qu.:    10.51  
>  Max.   : 999999    Max.   :103104.70  
>  NA's   :35820888  
>
>
> It seems that there are a lot of NA's. (Everything in the square outside the continent's borders should be an NA, so maybe that's a reasonable number). It seemed right when I plotted it earlier.
>
> I ran it with 9.333 as the distance band (then I was going to try 1.5 next).
> nb<- dnearneigh(r, 0, 9.333)
>
> It used a constant amount of memory for the first ~90 minutes and then it increased at a constant rate for the next 20 hours
> until it was using virtually all the memory on the server and was stopped.
> Because the administrator had to kill my job, I wasn't able to get session info while the session was open. I'm not sure this is at all helpful, but when I re-started R this is the sessionInfo() :
> R version 2.15.3 (2013-03-01)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C            LC_COLLATE=C        
>  [5] LC_MONETARY=C        LC_MESSAGES=C        LC_PAPER=C           LC_NAME=C           
>  [9] LC_ADDRESS=C         LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base     
>
> loaded via a namespace (and not attached):
> [1] tools_2.15.3
>
> Thanks,
> Juliann
>
>
> ________________________________
> From: Roger Bivand <Roger.Bivand at nhh.no>
> To: Empty Empty <phytophthorasb at yahoo.com> 
> Cc: "r-sig-geo at r-project.org" <r-sig-geo at r-project.org> 
> Sent: Sunday, April 7, 2013 12:15 PM
> Subject: Re: [R-sig-Geo] Memory problems with dnearneigh in spdep
> 
> I do prefer an affiliation, a real name, and a motivation; a
> yahoo or gmail address tells me nothing about what you might be expected
> to know. If I (or a list helper) cannot "place" you, it is much harder to
> know why you are having problems.
>
> You are also in breach of list rules because you are posting in HTML rather than just plain text. Removing the HTML will shorten your posting radically.
>
>
> On Sun, 7 Apr 2013, Empty Empty wrote:
>
>> Hi,
>> 
>> I'm running into memory problems calculating Gi* using dnearneigh and wonder if someone has any suggestions for either my script or another work around.
>> 
>> I have 1km data for Africa. My colleague calculated the appropriate distance band (9333m) using incremental spatial autocorrelation in ArcGIS but she was only able to run a subset of the data or with the data aggregated to 5km due to memory issues (fixed distance, Euclidean). I tried it in both ArcGIS (memory error after 45 hours) and R (Within 2 hours it was using more than 245GB on the server and then crashed). This is what I used
>
> Please always include the output of sessionInfo(). Here also include summary(r) output.
>
>> 
>> library(spdep)
>> 
>> 
>> r<-readOGR(".","r")
>> 
>> coords<-coordinates(r)
>
> This is not needed.
>
>> nb9333<- dnearneigh(coords, 0,9333,longlat = F)
>
> nb9333<- dnearneigh(r, 0, 9333)
>
> would work, and would not override longlat=. I assume in fact that r is in longlat, and that you would then mean to use a km threshold of 9.333. If I'm right, then the memory usage is being complicated by your including all observations as neighbours of each observation.
>
> For a 1km grid, a threshold of 9333m means that each of the about 30M gridcells has 100 neighbours (rather strong smoothing). Is the threshold sensibly chosen? The neighbour object will be 30M by 100 by 4, 1.21e+10 (12G); the labels of the neighbour object at least the same; the weights 24G in addition, so nb9333.listw is about 50G. With a more sensible threshold (1.5km?), the sizes of the nb9333.listw weights object become more tractable, maybe about 5G.
>
> Is any of the data observed on this grid? For localG* (or localG) to make any sense, the variable must be observed at this resolution (satellite imagery?) - if interpolated, all bets are off.
>
> Finally, the output involves so many multiple comparisons that no statistical inference is possible under normal assumptions - where is the 3.886 coming from - is this an attempt to adjust for multiple comparisons?
>
> I'm unsure why anyone would do this, as you could use standard filtering techniques on imaging data. The raster package also has functions for computing the local Moran and Geary measures on RasterLayer objects.
>
> Hope this clarifies,
>
> Roger
>
>> nb9333.listw <- nb2listw(include.self(nb9333), style="B")???
>> G9333<-localG(r$GRID_CODE,nb9333.listw)???
>> 
>> 
>> sig3 <- classIntervals(G9333, n=3, style="fixed",
>> ?????????????????????? fixedBreaks=c(min(G9333),-3.886, 3.886, max(G9333)))???? #I dont' think 3.886 is right
>> cs<-c('blue','white','red') # blue = below -3, white = intermediate, red = above 3
>> sigcol <- findColours(sig3, cs)?
>> 
>> png("r.png")
>> plot(r, col=sigcol, axes=TRUE)
>> dev.off()
>> 
>>     [[alternative HTML version deleted]]
>> 
>> 
>
> -- Roger Bivand
> Department of Economics, NHH Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand at nhh.no
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no