[R-sig-Geo] Maptools and readShapePoly
Roger Bivand
Roger.Bivand at nhh.no
Fri Sep 30 12:47:16 CEST 2005
On Fri, 30 Sep 2005, JP Azevedo wrote:
> Dear Roger,
>
> Many thanks for your attention on this matter.
>
> Last night I was trying to implement your suggestions but things did not get
> got as I expected. I've only manage to load the shape files using the old
> read.shape, however I did not manage to color the figure according to an
> auxvar. My priority is to try to use the readShapePoly.
>
> I'm running this from a Centrino 1.73 with 1.5 Gb of RAM.
>
Please also post the output of:
> version
> library(maptools)
> maptools()
> read.dcf(system.file("DESCRIPTION", package="sp"))[2]
so that we know what your system is.
>
> These are to examples of the readShapePoly syntax and both of them yield
> error messages.
>
> >
> >
> > brasilsp <- readShapePoly(system.file("shapes/Brasil_municipios_2001.shp",
> package="maptools")[1], proj4string=CRS("*"))
> Error in "dimnames<-.data.frame"(`*tmp*`, value = list(c("0", "1", "2", :
> invalid 'dimnames' given for data frame
> >
> > brasilsp <- readShapePoly("shapes\Brasil_municipios_2001.shp",
> proj4string=CRS("*"))
> Error in read.shape(filen = fn, verbose = verbose) :
> unable to open file
> >
I already suggested that you review where you keep your files. Please
never keep your own data inside the installed package directory tree, this
is quite unnecessary, and your files will be lost if you update the
package. Just keep your own shapefiles where you will be working (that is,
not where R is installed or where the packages are installed). On Windows,
you can easily change the working directory for R from the GUI File ->
menu, or use setwd(). On Windows you can also choose the shapefile
interactively by saying:
brasilsp <- readShapePoly(file.choose())
ignoring the proj4string unless you need it.
Once you can access the shapefile without putting it in harm's way in the
package own directory tree, we can get back to the
Error in "dimnames<-.data.frame"(`*tmp*`, value = list(c("0", "1", "2", :
invalid 'dimnames' given for data frame
question. If you can put the shapefile components in a zip archive on a
website, I can access them and see where that problem might be.
Roger
PS: the system.file() function is used in the examples to be able to
distribute shapefiles with the package, which may be installed by users in
varying locations. It is _not_ intended in any way to suggest that this
is the only place the function looks for data.
> >
> >
>
>
> When I use the old read.shape, I manage to load and plot the shape file if I
> use the full command, however the other two attempts fail.
>
> >
> > x <- read.shape(system.file("shapes/Brasil_municipios_2001.shp",
> package="maptools")[1])
> Shapefile type: Polygon, (5), # of Shapes: 5563
> >
> > x <- read.shape("shapes/Brasil_municipios_2001.shp", package="maptools")
> Error in read.shape("shapes/Brasil_municipios_2001.shp", package =
> "maptools") :
> unused argument(s) (package ...)
> >
> > x <- read.shape("shapes/Brasil_municipios_2001.shp")
> Error in read.shape("shapes/Brasil_municipios_2001.shp") :
> unable to open file
> >
> >
>
>
> Cheers,
>
> JP
>
> -----Original Message-----
> From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
> Sent: Friday, September 30, 2005 3:41 AM
> To: JP Azevedo
> Cc: R-sig-Geo at stat.math.ethz.ch
> Subject: Re: [R-sig-Geo] Maptools and readShapePoly
>
> 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
More information about the R-sig-Geo
mailing list