################################################### #### #### antipodes en projection orthogonale #### procrastin 2009 - procrastin@fxdm.org #### R version 2.8 #### #################################################### ##################### ## packages ##################### library(grDevices) library(maps) library(mapproj) library(mapdata) ##################### ## fonctions ##################### cercle <- function (xcentre, ycentre, rayon, lwd = 1, col = "black") { x <- seq(-rayon, rayon, by = 0.001) y <- sqrt(rayon^2 - x^2) lines(x+xcentre, y+ycentre, lwd = lwd, col = col) lines(x+xcentre, -y+ycentre, lwd = lwd, col = col) } circles <- function(x, y, r, col = rep(0, length(x)), border = rep(1, length(x)), ...) { for(i in 1:length(x)) { circle(x[i], y[i], r[i], col = col[i], border = border[i], ...) } } circle <- function(x, y, r, ...) { xx <- c() yy <- c() n <- 100 for(i in 1:n ) { ang <- i*6.28/n xx[i] <- cos(ang) yy[i] <- sin(ang) } xx <- r * xx + x yy <- r * yy + y polygon(xx, yy, ...) } mapfx <- function (database = "world", regions = ".", exact = FALSE, boundary = TRUE, interior = TRUE, projection = "", parameters = NULL, orientation = NULL, fill = FALSE, col = 1, plot = TRUE, add = FALSE, namesonly = FALSE, xlim = NULL, ylim = NULL, wrap = FALSE, resolution = if (plot) 1 else 0, type = "l", bg = par("bg"), mar = c(0, 0, par("mar")[3], 0.1), border = 0.01, mer = colmer, ...) { if (resolution > 0 && !plot) stop("must have plot=TRUE if resolution is given") if (!fill && !boundary && !interior) stop("one of boundary and interior must be TRUE") doproj <- !missing(projection) || !missing(parameters) || !missing(orientation) coordtype <- maptype(database) if (coordtype == "unknown") stop("missing database or unknown coordinate type") if (doproj && coordtype != "spherical") stop(paste(database, "database is not spherical; projections not allowed")) if (is.character(database)) as.polygon = fill else as.polygon = TRUE coord <- map.poly(database, regions, exact, xlim, ylim, boundary, interior, fill, as.polygon) if (is.na(coord$x[1])) stop("first coordinate is NA. bad map data?") if (plot) { assign(".map.range", coord$range, envir = globalenv()) } if (doproj) { nam <- coord$names library(mapproj) coord <- mapproject(coord, pr = projection, pa = parameters, or = orientation) coord$projection = projection coord$parameters = parameters coord$orientation = orientation if (plot && coord$error) if (all(is.na(coord$x))) stop("projection failed for all data") else warning("projection failed for some data") coord$names <- nam } if (plot) { if (!add) { opar = par(bg = bg) if (!par("new")) plot.new() if (is.null(xlim) || doproj) xrange <- range(coord$x, na.rm = TRUE) else xrange <- xlim if (is.null(ylim) || doproj) yrange <- range(coord$y, na.rm = TRUE) else yrange <- ylim if (coordtype != "spherical" || doproj) { aspect <- c(1, 1) } else aspect <- c(cos((mean(yrange) * pi)/180), 1) d <- c(diff(xrange), diff(yrange)) * (1 + 2 * border) * aspect if (coordtype != "spherical" || doproj) { plot.window(xrange, yrange, asp = 1/aspect[1]) circles(x= 0,y=0,r=1, col = mer) } else { p <- par("fin") - as.vector(matrix(c(0, 1, 1, 0, 0, 1, 1, 0), nrow = 2) %*% par("mai")) par(pin = p) p <- par("pin") p <- d * min(p/d) par(pin = p) d <- d * border + ((p/min(p/d) - d)/2)/aspect usr <- c(xrange, yrange) + rep(c(-1, 1), 2) * rep(d, c(2, 2)) par(usr = usr) } on.exit(par(opar)) } if (is.character(database) && resolution != 0 && type != "n") { pin <- par("pin") usr <- par("usr") resolution <- resolution * min(diff(usr)[-2]/pin/100) coord[c("x", "y")] <- mapthin(coord, resolution) } if (type != "n") { if (wrap) coord = map.wrap(coord) if (fill) polygon(coord, col = col, ...) else lines(coord, col = col, type = type, ...) } } class(coord) = "map" value <- if (namesonly) coord$names else coord if (plot) invisible(value) else value } ########################################## ## coordonnées du centre de projection ########################################## #### coordonnées de Paris centre_lat <- 48.51 centre_long <- 2.5 ####################### ### dessin du globe ####################### #### tracé du globe en vue orthogonale colmer = rgb(red = 153, green = 179, blue = 204, max = 255) colterre = rgb(red = 65, green = 108, blue = 108, max = 255) mapfx('world',proj='orth', orient=c(centre_lat, centre_long,0), resolution = 0, fill = T, col = colterre, mer = colmer) #### tracé des antipodes monde <- map('world', resolution = 0, fill = T, plot = F) antimonde <- monde antimonde$x <- monde$x + 180 antimonde$y <- -monde$y mapfx(antimonde,proj='orth',orient=c(centre_lat,centre_long,0), resolution = 0, col = rgb(red = 240, green = 240, blue = 240, max = 255, alpha = 100), fill = T, add= T) #### points rouge pour visualiser le centre de projection points(mapproject(list(y=centre_lat,x=centre_long)),col=2,pch=19,cex=0.5) #### dessin des méridiens, parallèles et tropiques colormeridien <- rgb(red = 110, green = 133, blue = 150, max = 255) nbr <- 6 # tracé des méridiens for (i in 1:nbr){ meridien <- seq(-90, 90, 0.01) x <- rep(2*180*i/nbr, length(meridien)) points(mapproject(list(y=meridien,x=x)),col=colormeridien, type = "l", lty = 1, lwd = 0.5) } # tracé des paralleles ecart = 40 latitude <- 0 while (latitude < 90) { parallele <- seq(-180, 180, 0.01) y <- rep(latitude, length(parallele)) points(mapproject(list(y=y,x=parallele)),col=colormeridien, type = "l", lty = 1, lwd = 0.5) latitude <- latitude + ecart } latitude <- 0 while (latitude > -90) { parallele <- seq(-180, 180, 0.01) y <- rep(latitude, length(parallele)) points(mapproject(list(y=y,x=parallele)),col=colormeridien, type = "l", lty = 1, lwd = 0.5) latitude <- latitude - ecart } # tracé des deux tropiques parallele <- seq(-180, 180, 0.01) y <- rep(23.5, length(parallele)) points(mapproject(list(y=y,x=parallele)),col=colormeridien, type = "l", lty = 3) y <- rep(-23.5, length(parallele)) points(mapproject(list(y=y,x=parallele)),col=colormeridien, type = "l", lty = 3) cercle(0,0,1, lwd = 3, col = "white")