[R-sig-Geo] Maptools and readShapePoly

[Roger Bivand]

On Thu, 29 Sep 2005, Roger Bivand wrote:

> On Thu, 29 Sep 2005, JP Azevedo wrote:
> 
> > Dear all,
> > 
> >  
> > 
> > I'm trying to use maptools to produce some figures using the Brazilian
> > municipalities shape files (5,500 municipalities and ~9mb).
> > 
> >  
> > 
> > I've ran some examples suggested by Roger which worked fine with the
> > sids.shp demo file.
> > 
> >  
> > 
> > library(maptools)
> > 
> >  
> > 
> > xsp <- readShapePoly(system.file("shapes/sids.shp",package="maptools")[1],
> > proj4string=CRS("+proj=longlat"))
> > 
> > plot(xsp)
> > 
> >  
> > 
> > brks <- seq(min(xsp$BIR74), max(xsp$BIR74), length=7) 
> > 
> > brks
> > 
> > cols <- colorRampPalette(c("lightpink", "red"))(length(brks)-1) 
> > 
> > cols
> > 
> >  
> > 
> > plot(xsp, col=cols[findInterval(xsp$BIR74, brks, all.inside=TRUE)])
> > 
> >  
> > 
> >  
> > 
> > However, whenever I try to load the Brazilian data, R just stalls, and
> > nothing happens.
> > 
> 
> "stalls" suggests that you have hit an error condition and that the 
> program has terminated abnormally. I think you just have to wait for it to 
> finish working. If you have any system monitor, you can track whether the 
> processor is fully occupied, and whether your disk is swapping (by ear).
> 
> >  
> > 
> > brasilsp <-
> > readShapePoly(system.file("shapes/Brasil_municipios_2001.shp",package="mapto
> > ols")[1], proj4string=CRS("+proj=longlat +ellps=clrk66"))
> > 
> 
> Well, using the very complicated system.file() function, which is only 
> needed for accessing, surprisingly, system files shipped with the maptools 
> package, makes things opaque.
> 
> > brasilsp <- readShapePoly("Brasil_municipios_2001.shp", 
> + proj4string=CRS("*")),
> 
> where * is the correct projection string for your shapefile, is easier. 
> Maybe the projection string isn't important for you now.
> 
> This will take a long time for 5000+ polygons, especially if the 
> boundaries are very detailed and have many islands and lakes. A map of 
> Australia with about 40K polygons and about 147 Mb was read successfully 
> (as far as I recall a longer lunch break on a regular office PC). I'll 
> have another try with the same shapefile on a 1.5GHz machine.

The test on the big Australian shapefile:

> version
         _                
platform i686-pc-linux-gnu
arch     i686             
os       linux-gnu        
system   i686, linux-gnu  
status                    
major    2                
minor    1.1              
year     2005             
month    06               
day      20               
language R       
> library(maptools)
Loading required package: foreign
Loading required package: sp
> system.time(bigshp <- readShapePoly("test1.shp"))
[1] 1952.79  118.77 2593.95    0.00    0.00
> object.size(bigshp)
[1] 216375656
> file.info("test1.shp")$size
[1] 147221660
> file.info("test1.dbf")$size
[1] 6940350
> np <- getSpPnParts(bigshp)
> length(np)
[1] 40349
> sum(np)
[1] 41209
> system.time(plot(bigshp))
[1]   94.77   25.71 1559.77    0.00    0.00

The R process ran to well over physical memory (1GB), so the difference 
between user time (the first number) and elapsed time (the third number) 
was swapping (my disk is not very fast). So the read took of the order av 
45 minutes, and could have been faster with more memory.

None of the small polygons were visible on the screen plot at all, with 
transparent borders, they might have got a pixel.

This is closer to your scenario:

> library(maptools)
Loading required package: foreign
Loading required package: sp
> system.time(bigshp <- readShapePoly("co99_d00.shp"))
[1] 131.30   0.61 134.28   0.00   0.00
> object.size(bigshp)
[1] 19088508
> file.info("co99_d00.dbf")$size
[1] 733012
> file.info("co99_d00.shp")$size
[1] 13108428
> np <- getSpPnParts(bigshp)
> length(np)
[1] 3489
> sum(np)
[1] 3525
> system.time(plot(bigshp))
[1] 4.53 0.32 5.43 0.00 0.00

which is for over 3000 counties in the US (continental and all islands, 
file from the US Census site). I think that is acceptable, given that the 
SpatialPolygonsDataFrame object can be saved as an R binary object and 
loaded when needed - it doesn't have to ve read from a shapefile every 
time.

Roger

> 
> However, very few useful statistical procedures can be done on this much 
> data, and visualization, which involved drawing each polygon separately, 
> is not going to be fast. R is not a GIS, and if visualising this kind of 
> data is your main concern, I would recommend Terraview, QGIS, or other 
> software aiming to get something to screen fast. 
> 
> Roger
> 
> >  
> > 
> > plot(brasilsp)
> > 
> >  
> > 
> > I've also tried to use the old  read.shape, which work fine the overall
> > plot.
> > 
> >  
> > 
> > x <- read.shape(system.file("shapes/Brasil_municipios_2001.shp",
> > package="maptools")[1])
> > 
> >  
> > 
> > plot(x, axes=FALSE, xlab = "", ylab = "")
> > 
> >  
> > 
> > plot(x, axes=FALSE, xlab = "", ylab = "", ol = "transparent")
> > 
> >  
> > 
> > However, if I try to use the auxvar option, I get an error message.
> > 
> >  
> > 
> > plot(x, auxvar=x$att.data$TOTAL, type="e",ol = "transparent", prbg=NULL)
> > 
> >  
> > 
> > Error in plot.Map(x, auxvar = x$att.data$TOTAL, type = "e", ol =
> > "transparent", :
> > 
> >  
> > 
> >  
> > 
> > Any help will be much appreciated.
> > 
> >  
> > 
> > Cheers,
> > 
> >  
> > 
> > JP
> > 
> >  
> > 
> > 
> > 	[[alternative HTML version deleted]]
> > 
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > 
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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