[R] adjusting the map of France to 1830
Stephane DRAY
dray at biomserv.univ-lyon1.fr
Fri Nov 19 01:24:21 CET 2004
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 : stephane.dray at umontreal.ca
--------------------------------------------------------------------------------------------------
Web http://www.steph280.freesurf.fr/
More information about the R-help
mailing list