adjusting the map of France to 1830

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/