Prevención del mapeo parcial de costas

Question

¿Cómo se puede evitar que los polígonos de un país se corten bajo diferentes proyecciones?

En el siguiente ejemplo, me gustaría hacer un mapa de proyección estereográfica de la Antártida, incluyendo las latitudes < -45 ° S. Al establecer mis límites y para este rango, el área de trazado es correcta, pero luego los polígonos de país también se recortan en estos límites. ¿Hay alguna manera de tener las costas trazadas en los bordes del área de trazado?

Gracias por cualquier consejo.

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

El problema se debe a que está especificando un máximo ylim de -45, y se están recortando los datos con precisión en esa latitud. Claramente, los límites reales de la trama se calculan de forma separada, lo que probablemente tenga que ver con algún tipo de suposición conservadora basada en los límites proyectados resultantes (pero de la que no sé nada sin explorar el código fuente).

Puede evitar el problema mediante la creación de una trama ficticia que no llama la línea de costa, a continuación, añadir los datos en la parte superior con un mayor 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() 

Por cierto, esta estrategia no funcionará para todas las proyecciones, ya que algunas implicarán una superposición para algunos límites de datos, pero Polar Stereographic simplemente se extiende desde el centro por lo que no hay tal problema aquí.

Además, habría una forma más actualizada de hacer esto con maptools, rgdal y rgeos que le permitirían utilizar el uso compartido de datos de una proyección PROJ.4 Stereographic adecuada (pero aún no he descifrado el recorte paso). Las proyecciones en mapproj son "solo de forma" y sería un poco difícil obtener otros datos en ellos, afaik.

Answer 2

Me di cuenta de que otra opción sería definir los límites de la región de trazado y no del mapa en sí. Esto podría dar cierta flexibilidad para definir el área de trazado exacta (es decir, xaxs = "i" y yaxs = "i"). También garantiza que se trazarán todos los polígonos: en el ejemplo de @mdsumner, Australia falta y los límites y deberían extenderse para trazar su polígono correctamente.

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

Otra forma es utilizar una proyección PROJ.4 real y utilizar el conjunto de datos en el paquete maptools. Aquí está el código, todavía hay un problema en el continente antártico donde la costa del océano se encuentra con el polo a -90 (eso es un poco molesto, y tomaría encontrar esa coordenada y eliminarla del resultado transformado, o ejecutar una herramienta de topología para encontrar la astilla).

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()