[R] adjusting the map of France to 1830

[Stephane DRAY]

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/