5

I have a shapefile of world countries, downloaded from here. I can plot it in R using

countries <- readOGR("shp","TM_WORLD_BORDERS-0.3",encoding="UTF-8",stringsAsFactors=F)
par(mar=c(0,0,0,0),bg=rgb(0.3,0.4,1))
plot(countries,col=rgb(1,0.8,0.4))

Now I wanna plot it in orthographic projection (Earth seen from outer space), so I'm trying

countries <- spTransform(countries,CRS("+proj=ortho +lat_0=-10 +lon_0=-60"))

I also played with the x_0 and y_0 parameters (as stated here), but I always get the error:

non finite transformation detected:
[1] 45.08332 39.76804      Inf      Inf
Erro em .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  : 
  failure in Polygons 3 Polygon 1 points 1
Além disso: Mensagens de aviso perdidas:
In .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  108 projected point(s) not finite

sometimes in the 3rd polygon, sometimes in the 7th. Where are those "Inf" coming from? I need to change any parameter? I want to plot the map like this

orthographic projection

but centered above South America. Thanks for your help!

Rodrigo
  • 4,706
  • 6
  • 51
  • 94
  • 1
    You are trying to project all the coordinates from -180 to 180 and from -90 to 90. But it is not possible with the kind of projection you are trying to use. Try first to crop to the region you want to see projected. And I strongly recommend you to read some literature about geographic projections. –  Sep 29 '14 at 14:25
  • See also: http://en.wikipedia.org/wiki/Orthographic_projection_in_cartography#Mathematics for a rule of thumb on where to crop. – plannapus Sep 29 '14 at 14:34
  • Thanks, Pascal. I wanna plot the whole half-world, just like the image above. I thought I would need only the central lat/long, and the algorithm would cut where needed. – Rodrigo Sep 29 '14 at 14:35
  • Thanks, plannapus, I'll try to use this formula to crop the points. Though I think that's just stupid. There should be an easier way to make this sort of plot automatically... – Rodrigo Sep 29 '14 at 14:38
  • In the first example (http://matplotlib.org/basemap/users/examples.html) they don't need to do what you're saying. It proves there's an easier way, or not? – Rodrigo Sep 29 '14 at 15:00
  • So @plannapus, you saw the link above? Don't you think that the function spTransform, since it calculates all the sin and cos, should already have an option to crop or not plot or anything else? Just like the example above seem to have? – Rodrigo Sep 30 '14 at 12:59
  • @Pascal? You saw the example above? – Rodrigo Sep 30 '14 at 12:59
  • Other point is, to crop the points hidden by Earth would destroy the polygons topology. – Rodrigo Oct 01 '14 at 14:19

1 Answers1

4

Give the maps package a try. It gives a warning about the points that cannot be projected, but it does not give an error and shut the process down. With a little fiddling, namely setting the fill color for the ocean (this answer helped solve that problem), I was able to emulate the map you attached with a couple lines:

library(maps)

## start plot & extract coordinates from orthographic map
o <- c(-10,-60,0) # oreantation
xy <- map("world",proj="orthographic", orientation=o, bg="black")
xy <- na.omit(data.frame(do.call(cbind, xy[c("x","y")])))

## draw a circle around the points for coloring the ocean 
polygon(max(xy$x)*sin(seq(0,2*pi,length.out=100)),max(xy$y)*cos(seq(0,2*pi,length.out=100)), 
        col="blue4", border=rgb(1,1,1,0.5), lwd=2)

## overlay world map
colRamp <- colorRampPalette(c("lemonchiffon", "orangered"))
map("world",proj="orthographic", orientation=o, 
    fill=TRUE, col=colRamp(5), add=TRUE)

enter image description here

Community
  • 1
  • 1
Paul Regular
  • 498
  • 4
  • 8
  • Thanks Paul, it is working. I'm doing it for an inset, giving the location of another, bigger map. How can I project the bounding box of the other map as a rectangle over this globe? (It will not be a rectangle, but a curved rectangle...) I could do it with mapproject and rect, but then I get a not curved rectangle. – Rodrigo Oct 02 '14 at 18:22
  • 1
    There are probably easier ways, but I would manually build the coords of a rectangle, project it using `mapproject` and add it using `polygon`: `x <- seq(-60, -30, length.out=100); y <- seq(-30, -10, length.out=100); x <- c(x, rep(x[length(x)], length(x)), rev(x), rep(x[1], length(x))); y <- c(rep(y[1], length(y)), y, rep(y[length(y)], length(y)), rev(y)); xy <- mapproject(x, y, proj="orthographic", orientation=o); polygon(xy$x, xy$y)` – Paul Regular Oct 02 '14 at 18:42
  • Yes, of course, how didn't I thought of that? ;) Thanks again! – Rodrigo Oct 02 '14 at 18:43