[R-sig-Geo] Maptools and readShapePoly

[Roger Bivand]

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.

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