在另一个国家内绘制许多国家 [英] Plot many countries inside another

问题描述

我想展示巴西亚马逊森林有多大，在里面绘制不同的国家.就像这张图片:

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

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

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

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

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/

按照@CiaPan 的回答，我来到了这个功能:

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)

}

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

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)

从此处下载 shapefile 的位置.谢谢大家！

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

推荐答案

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

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.

选择一个顶点，可能是最北端:V(φ0, λ0) 并决定它最终将落在亚马逊地区的何处:(φ1,λ1) 以及旋转多少:θ.您将通过几个简单的步骤实现它:

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. 沿着纬度圆滑动形状，使 V 落在格林威治子午线上——你可以从所有经度中减去 λ0:
   &λ;:= λ- λ0

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

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

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 φ

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

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)

将边框旋转 θ围绕 V 顶点，当前位于大西洋上的地理坐标 (0,0) – 这是一个平面 YZ 旋转:
Y := Y cos(θ) −Z sin(θ)
Z := Y sin(θ) + Z cos(θ)

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

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

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)

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

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

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

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

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

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

编辑

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

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.

这篇关于在另一个国家内绘制许多国家的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！