Skip to content
Prev 173718 / 398502 Next

Map using projection

OK, the underlying problem is that the projection code (i.e. the mapproj package plus some 
parts of the maps package) does not clip its output to the specified limits.

Let me explain:
1) the maps databases are made up of line segments which are combined to form polygons
2) when you call map() with limits specified, all line segements containing *at least one 
point* within the limits is returned
3) when there is no projection, the plot() function quite naturally clips the output to 
the specified rectangular limits
4) when there is a projection, the limits are no longer rectangular, and so this clipping 
no longer occurs, and all points are projected

Now with your code:
1) when you use the (default) "world" database, the line segments returned include a large 
part of the US-Canada border, plus all of Lake Winnipeg
2) if you use the "usa" database, you only get the US-Canada border but much more than 
that part specified by your limits
3)  if you use the "state" database, you get much closer to what you "want" (as defined by 
the limits), but you still get some parts of line segments which are outside the limits.

In order to get what you seem to want (although if you are looking at "states", why do you 
go such a long way north of the 49th parallel?), you can do something like this:

    library(maps)
    require("mapproj")
    longlatLimit <- c(-107, -93, 40, 52)
    par(plt=c(0, 1, 0, 1), cex=1, cex.main=1)  # Set plotting parameters
    mypoly <- list(x=c(-107, -93, -93, -107, -107), y=c(40, 40, 52, 52, 40))
    mymap <- map("state", plot=F, xlim=longlatLimit[1:2], ylim=longlatLimit[3:4])
    whichxy <- in.polygon(mypoly, list(x=mymap$x, y=mymap$y))
    whichxy[is.na(mymap$x)] <- NA
    mymap$x <- mymap$x[whichxy]
    mymap$y <- mymap$y[whichxy]
    map(mymap, projection="azequalarea", type="n")
    bound<-c(floor(longlatLimit[1]), ceiling(longlatLimit[2]),
      floor(longlatLimit[3]), ceiling(longlatLimit[4]))
    map.grid(lim=bound, col="light grey")
    map(mymap, projection="azequalarea", add=T)

HTH
Ray Brownrigg
Dr. Alireza Zolfaghari wrote: