At 19:07 18/11/2004, Ray Brownrigg wrote:
At 16:29 18/11/2004, Michael Friendly wrote: > I'm doing some analyses of historical data from France in 1830 on 'moral > statistics' that I'd like to > show on a map. I've done most of my analyses in SAS, but a few things > would work better in R. > To do this, I have to adjust the modern map, > > library(maps) > map('france') > > to adjust for changes in departments (86 in 1830, to 97 now). I've read > the documentation > for the maps and maptools package, but there seems to be no functions to > allow this, and > I can't find information on the exact structure of map datasets, but I > understand them to > be delimited lists of polygon coordinates. > > In SAS, all maps have (one or more) ID variables representing the > geographical region, > and there is also a proc gremove that can remove internal boundaries > inside the polygons > for regions with the same ID. Is there some way I can do this in R? > Unfortunately not with the current implementation of several of the 'extra' databases in the mapdata package. The map() function does have the interior=FALSE option, which would normally do what you want, but only when the data has been manipulated to allow it. Currently this option is only useful with the world and usa maps (and their derivatives, such as world2 and state).
Currently every department is a complete polygon, and so every interior line segment occurs twice in the data. What has to happen to the data is for it to be split up into non-overlapping line segments, and each polygon reconstructed from a list of these line segments (with direction being important).
There is the gpclib package which computes intersection, union... of polygons. I have try to play with its union function and the france data, but the results are good but a little bit complicate:
# i merge the two firts polygons:
> departements=map('france')
> which(is.na(departements$x))[1:2]
[1] 66 122
> gpcA <- as(cbind(departements$x[1:65],departements$y[1:65]),"gpc.poly")
> gpcB <- as(cbind(departements$x[67:121],departements$y[67:121]),"gpc.poly")
> union(gpcA,gpcB)
GPC Polygon
Num. Contours: 1
Num. Vertices: 74
BBox (X): 1.563161 --> 4.225965
BBox (Y): 49.97212 --> 51.09752
> gpcAB<-union(gpcA,gpcB)> departements$x=c(attr(gpcAB,"pts")[[1]]$x,attr(gpcAB,"pts")[[1]]$x[1],departements$x[-(1:121)])
> departements$y=c(attr(gpcAB,"pts")[[1]]$y,,attr(gpcAB,"pts")[[1]]$y[1],departements$y[-(1:121)])
> map(departements)
Another solution is to do the job in a GIS and to import the new map with maptools.
The package gpclib is very intersting and I think that it can be used to develop some basic GIS tools. I have write a functions to compute intersections for objects of class polys easily.
The only thing is to know for which class of spatial objects these functions must be developed !
Stéphane DRAY
--------------------------------------------------------------------------------------------------
Département des Sciences Biologiques Université de Montréal, C.P. 6128, succursale centre-ville Montréal, Québec H3C 3J7, Canada
Tel : (514) 343-6111 poste 1233 Fax : (514) 343-2293
E-mail : [EMAIL PROTECTED]
--------------------------------------------------------------------------------------------------
Web http://www.steph280.freesurf.fr/
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
