FW: [R] Drawing maps of UK

[Roger Bivand]

On Wed, 9 Feb 2005, michael watson (IAH-C) wrote:

> 
> Hi Guys
> 
> Thanks for this!  As I am a beginner, I bet I'm running into some really
> basic problems.  Using the example from the Map2poly function in
> maptools:
> 
> try <- read.shape("euadmll.shp")
> mappolys <-  Map2poly(try)
> # this produces 14 warnings() of the type:
> # 1: level 2 circularity at 62 in: .mtlbuild(.mtafters(r1), rD) # etc #
> 14: From next release, default hole handling will change in:
> Map2poly(try) 
> plot(mappolys)
> 
> Gives an error:
> Error in plot.polylist(mappolys) : malformed plot order
> 
> I also get this with Map2ploy(..., raw=FALSE)
> 
> I have no idea what this means! :-(

Basically, the original *.e00 file that Barry refered to had "topology", 
that is, polygons are built into rings on the fly from arcs, lines that 
are 2D coordinates. To draw a polygon, you need to build a ring getting 
the arcs in the right direction and order, otherwise you get spaghetti. 
Shapefiles just have rings, so can get pretty confused when the data they 
contain do not adhere to the ESRI standard (clockwise ring direction for 
islands, anti-clockwise for holes). Because most (unfortunately) 
shapefiles in the wild do not follow the specification, Map2poly() tries 
to work out whether ring directions seem plausible, and gives warnings 
when it gets confused. 

Finally. plot order is needed because R's polygon filling assumes that
polygons can be plotted one on top of the other, but this can lead to
wholly contained polygons (say Hamburg NUTS2 region)  being "painted over"
by its including polygon plotted later. Here the plot order heuristic
tries to see which plot order seems sensible, committing wholly contained
polygons last, or at least after the containing polygon if it is correctly
identified.

Downloading from: http://www.maths.lancs.ac.uk/~rowlings/Geog/Maps/

> library(maptools)
Loading required package: foreign 
> try <- read.shape("euadmll.shp")
Shapefile Type: Polygon   # of Shapes: 2556 
> mappolys <-  Map2poly(try, raw=FALSE)
There were 50 or more warnings (use warnings() to see the first 50)

One possible problem with the file is that the first shape seems to be the 
whole of Europe as a single polygon:

>  str(mappolys[[1]])
 num [1:202333, 1:2] 25.6 25.5 25.5 25.5 25.4 ...
 - attr(*, "nParts")= int 2003
 - attr(*, "pstart")=List of 2
  ..$ from: int [1:2003] 1 130 185 198 260 539 28921 29056 29138 29159 ...
  ..$ to  : int [1:2003] 128 183 196 258 537 28919 29054 29136 29157 29173 
...
 - attr(*, "bbox")= num [1:4] -24.5  27.6  31.6  71.2
 - attr(*, "ringDir")= num [1:2003] 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "after")= int [1:2003] NA NA NA NA NA NA NA NA NA NA ...
 - attr(*, "plotOrder")= int [1:2003] 1 2 3 4 5 6 190 7 8 9 ...
 - attr(*, "shpID")= int NA

and this seems to be the problem. This does happen in e00 files - the
first object is the totality of the objects included, giving an outer
border, here with 2003 islands, lakes, etc. The plot order malformation is
then clear, everything is contained inside itself in the first shape. A
possible fix is then:

> length(mappolys)
[1] 2556
> mp1 <- subset(mappolys, 1:length(mappolys) %in% 2:length(mappolys))
> plot(mp1)

which works for me.

> names(try$att.data)
[1] "AREA"      "PERIMETER" "NUM"       "ID"        "COUNTRY"   "TOTALAREA"

suggests though that identifying which bits of the jigsaw-puzzle are which 
may not be easy, perhaps there is more to go on in the original e00.

Roger

> 
> Mick
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no