将许多国家列为另一个国家 [英] Plot many countries inside another
问题描述
为了实现这一点，我加载了一些shapefile，并将其投影更改为将保持区域成比例的投影，如圆柱平均面积：
library（rgdal）
国家< - readOGR（shp，TM_WORLD_BORDERS-0.3）
国家< - spTransform ，CRS（+ proj = cea））
amzLegal< - readOGR（shp，amazlegal）
amzLegal @ proj4string< - CRS（+ proj = longlat $ b amzLegal< - spTransform（amzLegal，CRS（+ proj = cea））
plot（amzLegal）
FR< - 国家[哪些（国家$ NAME ==法国）， ]
for（i in 1：length（FR @ polygons [[1]] @ Polygons））{
FR @ polygons [[1]] @ Polygons [[i]] @ coords [，1 ] = FR @多边形[[1]] @多边形[[i]] @ coords [，1] -7180000
FR @多边形[[1]] @多边形[[i]] @ coords [，2] = FR @多边形[[1]] @多边形[[i]] @ coords [，2] -4930000
}
plot（FR，col =blue，add = T）
我得到这个（没有线，我稍后添加）：
根据Google地球，红线约950公里（在法国），黑线一样（巴西）。所以当然圆柱形平均面积不是正确的投影，因为它扩大了经度，缩小了纬度。那么我应该用什么投影呢？一个保持形状和大小？我也尝试了兰伯特方位角平衡区，但也没有工作。我喜欢古德的Homolosine，但它不是一个单一的投影，而是混合不同的技术。以下是可能的预测列表： http://www.remotesensing.org/geotiff/proj_list/
编辑：以下@CiaPan答案，我来到这个功能：
translate < - function（obj，x，y，ang = 0，adiciona = T）{
maxLat <-90
for（i in 1 ：long（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＃度到弧度
ang< - ang * pi / 180
x1 = 180
x2 = -180
y1 = 90
y2 = -90
for（i in 1：length（obj @ polygons [ 1]] @多边形））{
for（j in 1：nr ow（obj @ polygons [[1]] @ Polygons [[i]] @ coords））{
lon< - obj @ polygons [[1]] @ Polygons [[i]] @ coords [ 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笛卡尔坐标
Y < - sin（lon）* cos（lat）
Z < ）
X0 < - X
X <-X0 * cos（lat0） - Z * sin（-lat0）＃3 V到赤道
Z <-X0 * sin（-lat0）+ Z * cos（lat0）
Y0 < - Y
Y < - Y0 * cos（ang）-Z * sin（ang） 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}
如果（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） ）
国家< - readOGR（shp，TM_WORLD_BORDERS-0.3，encoding =UTF-8）
plot（countries，col = rgb（1,0.8,0.4））
translate（国家（国家$ NAME ==法国），]， - 60,0,0，T）
其中shapefile是从 here 下载的。谢谢大家！
首先假设你的国家的边界被赋予地理（&phi，λ）坐标 - 如果他们在某些制图投影中是（x，y），你必须将它们转换回地理系统。
选择一个顶点，可能是最北端的：V（&phi ; 0，λ 0），并决定它将最终落在亚马逊地区：（φ 1，λ 1）和旋转多少：θ你可以通过几个简单的步骤来实现：
-
沿着纬度的圆圈滑动形状，使V落在格林威治子午线 - 你从所有经度中减去λ 0;
λ ：=λ - λ 0 -
接下来计算滑动边界的所有顶点的笛卡尔坐标（假设地球表面是球体，而不是椭圆体，并将地球半径作为长度单位）：
X：= cosλ cosφ
Y：= sinλ cosφ
Z：= sinφ -
将形状向南滑动，使V落在赤道上。您可以通过（−φ 0）角度旋转XZ平面中的所有顶点：
X：= X cos（φ 0）− Z sin（−φ 0）
Z：= X sin（−φ 0）+ Z cos（φ 0） -
通过θ旋转边框在当前位于大西洋的地理坐标（0,0）的V顶点附近，这是一个平面YZ旋转：
Y：= Y cos（θ）− Z sin（θ）
Z：= Y sin（θ）+ Z cos（θ） -
边界准备去亚马逊森林旅行。首先沿着格林威治子午线向南滑动到所需的纬度（平面XZ旋转φ 1 - noteφ 1为负，因为它表示南半球）：
X：= X cos（&phi ; 1）− Z sin（φ 1）
Z：= X sin（φ 1）+ Z cos（φ 1） -
然后将坐标转换为地理系统：
φ ：= asin（Z）
λ ：= asin（Y / cos（&phi）） -
最后将它们滑到南美洲
λ =λ +λ 1 -
完成。至少我希望如此...;）
编辑
您还可以在7.之后的1.和6.之前执行步骤2.
然后，当然，沿着纬度的圆滑动边框不会和lambda一样简单：=λ + const，它必须被计算为XY平面旋转，类似于步骤3到5.然而，这样一来，所有的变换将以类似的方式执行，您可以将其描述为矩阵乘法。矩阵乘法是相关的，所以所有系数矩阵都可以预先计算并乘以（按正确顺序！），然后用单个矩阵乘法转换边界的每个顶点。
一旦你处理了所有的国家，只是绘制所有的国家，看看它们是否相交。在这种情况下，调整目标点和θ;旋转，直到所有的边界适合亚马逊的丛林轮廓，没有碰撞。希望有帮助。
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!
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:
Slide the shape along the circle of latitude, so that V lands on the Greenwich meridian – you do this subtracting λ0 from all longitudes:
λ := λ - λ0Next 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 φ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)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(θ)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)then transform coordinates to geographic system:
φ := asin(Z)
λ := asin(Y/cos(φ))and finally slide them west to the South America
λ = λ + λ1Done. 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.
这篇关于将许多国家列为另一个国家的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
