[R] Map using projection

[Ray Brownrigg]

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:
> Ray,
> You are right, it is not having any data for states, but in fact the 
> grids are related to states. If you run the code you sent, you could see 
> there is a blank area on top and left. All I want is to locate the grid 
> map in the middle of window. Once I do this I can add the boundary data.
>  
> Thanks,
> Alireza
> 
> On Fri, Mar 13, 2009 at 9:00 AM, Ray Brownrigg 
> <Ray.Brownrigg at ecs.vuw.ac.nz <mailto:Ray.Brownrigg at ecs.vuw.ac.nz>> wrote:
> 
>     I think what you are missing is that the default map database is
>     "world".
> 
>     If you use:
> 
>        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
>        map(regions="USA", projection="azequalarea",
> 
>     type="n",xlim=longlatLimit[1:2],ylim=longlatLimit[3:4])
>        bound<-c(floor(longlatLimit[1]), ceiling(longlatLimit[2]),
>     floor(longlatLimit[3]), ceiling(longlatLimit[4]))
>        map.grid(lim=bound,col="light grey")
> 
>     I think you will get what you want (but what you want is not exactly
>     clear).
> 
>     [I am not sure why this happens, but it seems to be an artifact of
>     the projection code.]
> 
>     [All your examples have been very confusing, since you talk about
>     plotting the "states", but in fact your example code plots only a grid.]
> 
>     HTH
>     Ray Brownrigg
> 
> 
>     Dr. Alireza Zolfaghari wrote:
> 
>         Hi list,
>         I have a real problem with plotting US state map. When I try to
>         plot the
>         northern state, there will be some blank space in the top of
>         graph (see case
>         1 example), and when I plot southern states, there will be a
>         blank space in
>         the bottom of plot (see case 2). I spent almost 2 days to figure
>         out a
>         solution, but could not. Would you help me if you know what the
>         problem is?
>         Regards,
>         Alireza
> 
>         #################################################
>         #case 1
>            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
>            #map(projection="azequalarea",
>         type="n",xlim=longlatLimit[1:2],ylim=longlatLimit[3:4])
>            map(projection="azequalarea",
>         type="n",xlim=longlatLimit[1:2],ylim=longlatLimit[3:4])
>            bound<-c(floor(longlatLimit[1]), ceiling(longlatLimit[2]),
>         floor(longlatLimit[3]), ceiling(longlatLimit[4]))
>            map.grid(lim=bound,col="light grey")
>         #################################################
>         #case 2
>            library(maps)
>            require("mapproj")
>            longlatLimit<-c(-107,-93,25,37)
>            par(plt=c(0,1,0,1),cex=1,cex.main=1)  #Set plotting parameters
>            #map(projection="azequalarea",
>         type="n",xlim=longlatLimit[1:2],ylim=longlatLimit[3:4])
>            map(projection="azequalarea",
>         type="n",xlim=longlatLimit[1:2],ylim=longlatLimit[3:4])
>            bound<-c(floor(longlatLimit[1]), ceiling(longlatLimit[2]),
>         floor(longlatLimit[3]), ceiling(longlatLimit[4]))
>            map.grid(lim=bound,col="light grey")
> 
>                [[alternative HTML version deleted]]
> 
>         ______________________________________________
>         R-help at r-project.org <mailto:R-help at r-project.org> mailing list
>         https://stat.ethz.ch/mailman/listinfo/r-help
>         PLEASE do read the posting guide
>         http://www.R-project.org/posting-guide.html
>         <http://www.r-project.org/posting-guide.html>
>         and provide commented, minimal, self-contained, reproducible code.
> 
> 
> 
>