在另一个内部绘制许多国家

Question

我想展示巴西亚马逊森林有多大，在其中绘制不同的国家。就像这张图片：

为了实现这一点，我加载了一些 shapefile 并将它们的投影更改为保持区域成比例的投影，例如圆柱等面积：

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)

我得到了这个（不包括我稍后添加的行）：

根据 Google 地球，红线约为 950 公里（在法国），与黑线（在巴西）的长度相同。因此，圆柱等面积当然不是合适的投影，因为它扩大了经度并缩小了纬度。那我应该使用什么投影呢？一个保持形状和大小的？我也尝试过 Lambert Azimuthal Equal Area，但也没有用。我喜欢 Goode 的 Homolosine，但它并不是真正的单一投影，而是不同技术的混合。以下是可能的预测列表：http://www.remotesensing.org/geotiff/proj_list/

编辑：在@CiaPan 回答之后，我来到了这个函数：

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)

}

其中 obj 是一个 spatialPolygons 对象，x 和 y 是目标的 long 和 lat。该函数翻译并绘制对象。用法可以是：

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)

shapefile 是从here 下载的。谢谢大家！

Answer 1

首先假设您的国家的边界​​是用地理 (φ,λ) 坐标给出的 - 如果它们在某些制图投影中是 (x,y)，您必须将它们转换回地理系统。

选择一个顶点，可能是最北端：V(φ0, λ0) 并决定它最终在亚马逊地区的位置：(φ1,λ1) 和旋转多少：θ。您将通过几个简单的步骤来实现它：

1. 沿纬度圈滑动形状，使 V 落在格林威治子午线上 - 你这样做是从所有经度中减去 λ0：
   λ := λ - λ0

2. 接下来计算滑动边框所有顶点的笛卡尔坐标（假设地球表面是球体，不是椭球更不用说大地水准面，以地球半径为长度单位）：
   X := cos λ cos φ
   Y := sin λ cos φ
   Z := sin φ

3. 将图形向南滑动，使 V 落在赤道上。您可以将 XZ 平面中的所有顶点旋转 (−φ0) 角：
   X := X cos(φ0) − Z sin(−φ0)
   Z := X sin(−φ0) + Z cos(φ0)

4. 围绕 V 顶点将边界旋转 θ，该顶点当前位于大西洋的地理坐标 (0,0) - 这是平面 YZ 旋转：
   Y := Y cos(θ) - Z sin(θ)
   Z := Y sin(θ) + Z cos(θ)

5. 现在国家边界已准备好前往亚马逊森林。首先沿着格林威治子午线向南滑动到所需的纬度（平面 XZ 旋转 φ1 - 注意 φ1 为负，因为它表示南半球）：
   X := X cos(φ1) − Z sin(φ1)
   Z := X sin(φ1) + Z cos(φ1)

6. 然后将坐标转换为地理系统：
   φ := asin(Z)
   λ := asin(Y/cos(φ))

7. 最后将它们向西滑到南美洲
   λ = λ + λ1

8. 完成。至少我希望如此... ;)

编辑

您也可以在 1 之前执行第 2 步，在 7 之后执行第 6 步。
然后，当然，沿纬度圆滑动边界不会像 λ := λ + const 那样简单，它必须计算为 XY 平面旋转，类似于步骤 3 到 5。这样，但是，所有转换都将以类似的方式执行，您可以将其描述为矩阵乘法。并且矩阵乘法是关联的，因此可以预先计算所有系数矩阵并相乘（以正确的顺序！），然后用单个矩阵乘法转换边界的每个顶点。

处理完所有国家/地区后，只需将它们全部绘制出来，看看它们是否相交。在这种情况下，调整目标点和 θ 旋转，直到所有边界都符合亚马逊丛林轮廓而没有碰撞。希望对您有所帮助。

Answer 2

我建议尝试使用不同的投影，然后在您生成的地图上绘制 Tissot 的椭圆（下面的前两个链接）。您可以直观地检查地图并选择具有相似失真的国家。

如果您只是想在视觉上进行比较，任何中断的投影都是最好的。唯一的问题是有很多不连续性。每次你想创建一个国家的图像时，你都会改变投影，直到整个国家（或尽可能多的）没有间断。 仅通过浏览您的列表，我没有看到任何我认为被打断的内容。如果您不受这些预测的严格限制，我推荐 Goode's Homolosine，因为它将不连续性置于海洋中。

参考：

• 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

此软件（免费）允许您在许多不同的投影上进行比较（并绘制 tissot 的椭圆）： http://www.flexprojector.com/