Prevent partial display of coastlines

Question

How to prevent disconnection of country polygons in different projections?

In the following example, I would like to make a stereographic projection map of Antarctica, including latitudes <-45 ° C. Having set my y limits for this range, the graph area is correct, but then the country's polygons are also cut off within these limits. Is there a way for coastlines to be built around the edges of the plot area?

Thanks for the consultation.

library(maps) library(mapproj) ylim=c(-90,-45) orientation=c(-90, 0, 0) x11() par(mar=c(1,1,1,1)) m <- map("world", plot=FALSE) map("world",project="stereographic", orientation=orientation, ylim=ylim) map.grid(m, nx=18,ny=18, col=8) box() 

Answer 1

The problem is that you specify max ylim of -45, and the data is precisely cut at that latitude. Obviously, the actual boundaries of the graph are calculated in a separate way, which is probably due to some conservative assumption based on the obtained predicted boundaries (but I don’t know anything about it without studying the source code).

You can avoid the problem by setting a dummy chart that does not draw the coastline, then add the data at the top with a larger ylim :

 library(maps) library(mapproj) ylim = c(-90,-45) orientation=c(-90, 0, 0) x11() par(mar=c(1,1,1,1)) m <- map("world", plot=FALSE) map("world",project="stereographic", orientation=orientation, ylim=ylim, col = "transparent") map("world",project="stereographic", orientation=orientation, ylim=ylim + c(-5, 5), add = TRUE) map.grid(m, nx=18,ny=18, col=8) box() 

By the way, this strategy will not work for all forecasts, since some of them assume overlapping for some data limits, but Polar Stereographic just goes out of the center, so there is no such problem.

In addition, there would be a more modern way to do this with maptools, rgdal and rgeos, which will allow you to use data to share PROJ.4 Stereographic projection (but I did not quite understand the clipping step). Projections in mapproj are “only formal”, and it would work a bit to get other data in afaik ones.

Answer 2

I realized that another option would be to determine the boundaries of the chart area, not the map itself. This may provide some flexibility in determining the exact area of ​​the graph (ie, Xaxs = "i" and yaxs = "i"). You also guarantee that all polygons will be built - in the example @mdsumner, Australia is absent, and the y-limits will need to be expanded in order to correctly build the polygon.

 orientation=c(-90, 0, 0) ylim <- c(mapproject(x=-180,y=-45, project="stereographic", orientation=orientation)$y, mapproject(x=0,y=-45, project="stereographic", orientation=orientation)$y) xlim <- c(mapproject(x=-90,y=-45, project="stereographic", orientation=orientation)$x, mapproject(x=90,y=-45, project="stereographic", orientation=orientation)$x) x11(width=6,height=6) par(mar=c(1,1,1,1)) plot(0,0, t="n", ylim=ylim, xlim=xlim, xaxs="i", yaxs="i", xlab="", ylab="", xaxt="n", yaxt="n") map("world",project="stereographic", orientation=orientation, add=TRUE) map.grid(nx=18,ny=18, col=8) box() 

Answer 3

Another way is to use the actual PROJ.4 projection and use the dataset in the maptools package. Here's the code, there is still a problem on the Antarctic continent where the coastline of the ocean meets the -90 pole (which is a little annoying, and he will find this coordinate and remove it from the converted result or execute a topology tool to find the tape).

 library(maptools) data(wrld_simpl) xrange <- c(-180, 180, 0, 0) yrange <- c(-90, -45, -90, -45) stere <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" library(rgdal) w <- spTransform(wrld_simpl, CRS(stere)) xy <- project(cbind(xrange, yrange), stere) plot(w, xlim = range(xy[,1]), ylim = range(xy[,2])) box()