Plot many countries inside another

问题:

I want to show how big Brazilian Amazon Forest is, plotting different countries inside it. Like in this image:

To accomplish that, I loaded some shapefiles and changed their projection to one that would keep the areas proportional, like Cylindrical Equal Area:

library(rgdal)
countries <- readOGR("shp","TM_WORLD_BORDERS-0.3")
countries <- spTransform(countries,CRS("+proj=cea"))
amzLegal <- readOGR("shp","amazlegal")
amzLegal@proj4string <- CRS("+proj=longlat")
amzLegal <- spTransform(amzLegal,CRS("+proj=cea"))
plot(amzLegal)
FR <- countries[which(countries$NAME == "France"),]
for (i in 1:length(FR@polygons[[1]]@Polygons)) {
FR@polygons[[1]]@Polygons[[i]]@coords[,1] = FR@polygons[[1]]@Polygons[[i]]@coords[,1]-7180000
FR@polygons[[1]]@Polygons[[i]]@coords[,2] = FR@polygons[[1]]@Polygons[[i]]@coords[,2]-4930000
}
plot(FR,col="blue",add=T)

I'm getting this (without the lines, which I added later):

According to Google Earth, the red line is about 950 km (in France), the same measure of the black line (in Brazil). So of course the Cylindrical Equal Area is not the proper projection to use, since it enlarges the longitude and shrinks the latitude. What projection should I use, then? One that keeps shape AND size? I have also tried the Lambert Azimuthal Equal Area, but didn't work either. I like the Goode's Homolosine but it's not really a single projection, but a mix of different techniques. Here is a list of the possible projections: http://www.remotesensing.org/geotiff/proj_list/

EDIT: Following @CiaPan answer, I came to this function:

translate <- function(obj,x,y,ang=0,adiciona=T) {

maxLat <- -90
for (i in 1:length(obj@polygons[[1]]@Polygons)) {
    for (j in 1:nrow(obj@polygons[[1]]@Polygons[[i]]@coords)) {
        lat <- obj@polygons[[1]]@Polygons[[i]]@coords[j,2]
        if (lat > maxLat) {
            maxLat <- lat
            maxLon <- obj@polygons[[1]]@Polygons[[i]]@coords[j,1]
        }
    }
}
lon0 <- maxLon*pi/180
lat0 <- maxLat*pi/180

y <- y*pi/180 # degrees to radians
ang <- ang*pi/180
x1 = 180
x2 = -180
y1 = 90
y2 = -90
for (i in 1:length(obj@polygons[[1]]@Polygons)) {
    for (j in 1:nrow(obj@polygons[[1]]@Polygons[[i]]@coords)) {
        lon <- obj@polygons[[1]]@Polygons[[i]]@coords[j,1]*pi/180 - lon0      #1 V to Greenwich
        lat <- obj@polygons[[1]]@Polygons[[i]]@coords[j,2]*pi/180

        X <- cos(lon)*cos(lat)                 #2 Cartesian coords
        Y <- sin(lon)*cos(lat)
        Z <- sin(lat)

        X0 <- X
        X <- X0*cos(lat0) - Z*sin(-lat0)       #3 V to Equator
        Z <- X0*sin(-lat0) + Z*cos(lat0)

        Y0 <- Y
        Y <- Y0*cos(ang) - Z*sin(ang)          #4 rotate by ang
        Z <- Y0*sin(ang) + Z*cos(ang)

        X0 <- X
        X <- X0*cos(y) - Z*sin(y)              #5 V to y
        Z <- X0*sin(y) + Z*cos(y)

        lat <- asin(Z)                         #6
        lon <- asin(Y/cos(lat))*180/pi + x
        lat <- lat*180/pi

        if (lon < x1) { x1 <- lon }            #bbox
        if (lon > x2) { x2 <- lon }
        if (lat < y1) { y1 <- lat }
        if (lat > y2) { y2 <- lat }

        obj@polygons[[1]]@Polygons[[i]]@coords[j,1] <- lon
        obj@polygons[[1]]@Polygons[[i]]@coords[j,2] <- lat
    }
}
obj@bbox[1,1] <- x1
obj@bbox[1,2] <- x2
obj@bbox[2,1] <- y1
obj@bbox[2,2] <- y2

plot(obj,col="red",border="black",add=adiciona)

}

Where obj is a spatialPolygons object, x and y are long and lat of destination. The function translates and plot the object. Usage can be:

library(rgdal)
par(mar=c(0,0,0,0))
countries <- readOGR("shp","TM_WORLD_BORDERS-0.3",encoding="UTF-8")
plot(countries,col=rgb(1,0.8,0.4))
translate(countries[which(countries$NAME == "France"),],-60,0,0,T)

where the shapefile was downloaded from here. Thank you all!

回答1:

First assuming your countries' borders are given with geographic (φ,λ) coordinates – if they are (x,y) in some cartographic projection, you'd have to transform them back to the geographic system.

Choose one vertex, possibly the northernmost: V(φ0, λ0) and decide where it shall finally land in the amazonian region: (φ1,λ1) and how much rotated: θ. You'll achieve it in several simple steps:

1. Slide the shape along the circle of latitude, so that V lands on the Greenwich meridian – you do this subtracting λ0 from all longitudes:
   λ := λ - λ0

2. Next calculate Cartesian coordinates of all vertices of the slided border (assuming the Earth surface is a sphere, not an ellipsoid let alone the geoid, and taking the Earth radius as a length unit):
   X := cos λ cos φ
   Y := sin λ cos φ
   Z := sin φ

3. Slide the shape to the south, so that V lands on the equator. You do this rotating all vertices in XZ plane by the (−φ0) angle:
   X := X cos(φ0) − Z sin(−φ0)
   Z := X sin(−φ0) + Z cos(φ0)

4. Rotate the border by θ around the V vertex, which currently lays on Atlantic Ocean at geographic coordinates (0,0) – this is a plane YZ rotation:
   Y := Y cos(θ) − Z sin(θ)
   Z := Y sin(θ) + Z cos(θ)

5. Now the country border is ready to travel to Amazonian Forests. First slide it south along the Greenwich meridian to the desired latitude (plane XZ rotation by φ1 – note φ1 is negative, as it denotes the southern hemisphere):
   X := X cos(φ1) − Z sin(φ1)
   Z := X sin(φ1) + Z cos(φ1)

6. then transform coordinates to geographic system:
   φ := asin(Z)
   λ := asin(Y/cos(φ))

7. and finally slide them west to the South America
   λ = λ + λ1

8. Done. At least I hope so... ;)

EDIT

You also might do step 2. before 1. and step 6. after 7.
Then, of course, sliding the border along the circle of latitude would not be as simple as λ := λ + const., it would have to be computed as a XY-plane rotation, similar to steps 3. through 5. This way, however, ALL transformations will be performed in a similar way, which you can describe as a matrix multiplication. And matrix multiplication is associative, so all coefficient matrices can be calculated in advance and multiplied together (in the proper order!), then you transform each vertex of a border with a single matrix multiplication.

Once you processed all the countries just plot them all to see if they intersect. In such case tune the destination points and θ rotations until all borders fit the Amazonian jungle contour with no collision. Hope that helps.

回答2:

I recommend using trying different projections, but then plot Tissot’s ellipses on the maps you generate (first two links below). You could visually inspect the maps and pick countries that have similar distortion.

If you just want to compare visually, any interrupted projection would be best. The only problem is that there are a lot of discontinuities. Each time you want to create an image of a country, you shift the projection until the entire country (or as much as you can get) is without discontinuity. Just from scanning through your list, I don't see any that I recognize as interrupted. If you are not strictly limited to these projections, I recommend Goode's Homolosine as it puts the discontinuities in oceans.

Ref:

• http://en.wikipedia.org/wiki/Tissot%27s_Indicatrix
• https://gis.stackexchange.com/questions/5068/how-to-create-an-accurate-tissot-indicatrix
• http://proceedings.esri.com/library/userconf/proc98/PROCEED/TO850/PAP844/P844.HTM

This software (free), allows you to compare (& plot tissot's ellipses) on many different projections: http://www.flexprojector.com/