在绘制世界地图时使用与主要子午线不同的中心 [英] Use different center than the prime meridian in plotting a world map
问题描述
我将 maps
包中的世界地图覆盖到 ggplot2
栅格几何上。然而,这个栅格不是以主子午线(0度)为中心,而是在180度(大约白令海和太平洋)。以下代码获取地图并重新映射180度地图:
require(地图)
world_map = data。 frame(map(plot = FALSE)[c(x,y)])
names(world_map)= c(lon,lat)
world_map = within(world_map, {
lon = ifelse(lon <0,lon + 360,lon)
})
ggplot(aes(x = lon,y = lat),data = world_map)+ geom_path( )
产生以下输出:
直线在主子午线一端或另一端的多边形之间绘制。我目前的解决方案是用NA代替接近主子午线的点,用中的调用替换:
<$ p (世界地图,{
lon = ifelse(lon <0,lon + 360,lon)
lon = ifelse((lon <1))内的世界地图。 (lon> 359),NA,lon)
})
ggplot(aes(x = lon,y = lat),data = world_map)+ geom_path()
这导致了正确的图像。我现在有许多问题:
- 在另一条子午线上定位地图必须有更好的方法。我尝试在
map
中使用orientation 参数,但将其设置为
orientation = c(0,180 ,0)
没有产生正确的结果,事实上它没有改变任何东西到结果对象(all.equal
产生TRUE
)。 - 在不删除某些多边形的情况下,应该可以摆脱水平条纹。这可能是解决方法1。也解决了这一问题。 。它的工作原理是:
- 将世界地图从
maps
包转换为SpatialLines
具有地理(lat-long)CRS的对象。 - 投影
SpatialLines
映射到以素数子午线为中心的PlateCarée(aka等距圆柱形)投影。 (此投影与地理映射非常相似)。
- 剪切两段,否则将被地图的左右边缘剪切。 (这是使用来自
rgeos
包的拓扑函数完成的。) - 重新投影到PlateCarée投影,集中于所需子午线(<从
PROJ_4
程序中使用的术语中,code> lon_0spTransform()
在rgdal
包中)。 - 识别(并删除)任何剩余的条纹。我通过搜索通过g.e.的线来自动执行此操作。三个分开的经络中的两个。 (这也使用来自
rgeos
包的拓扑函数。)
这是显然有很多工作要做,但留下一个最小截断的地图,并且可以使用
spTransform()
轻松重新投影。要使用base
或lattice
图形将这些覆盖在光栅图像上,我首先重新映射栅格,同时使用spTransform()
。如果你需要它们,网格线和标签也可以预计为匹配SpatialLines
地图。
库(sp)
库(映射)
库(maptools)## map2SpatialLines(),pruneMap()
库(rgdal)## CRS(),spTransform()$ b $库库(rgeos)## readWKT(),gIntersects(),gBuffer(),gDifference()
## map映射到SpatialLines映射
makeSLmap< - function(){
llCRS< - CRS(+ proj = longlat + ellps = WGS84)
wrld< - map(world,interior = FALSE,plot = FALSE,
xlim = c(-179,179),ylim = c(-89,89))
wrld_p < - pruneMap(wrld ,xlim = c(-179,179))
map2SpatialLines(wrld_p,proj4string = llCRS)
}
## Clip SpatialLines整齐地沿着反对子午线
sliceAtAntipodes < - 函数(SLmap,lon_0){
##预备
long_180< - (lon_0 %% 360) - 180
llCRS< - CRS(+ proj = longlat + ellps = 地图对象的## CRS
eqcCRS < - CRS(+ proj = eqc)
##将地图重新映射到Equidistant圆柱/板块Caree投影
SLmap< ; - 线(cbind(long_180,c(-89,89))),ID(线)(cbind(long_180,c(-89,89))); - 在子午线上形成一个窄的SpatialPolygon =cutter)
SL < - SpatialLines(list(L),proj4string = llCRS)
SP < - gBuffer(spTransform(SL,eqcCRS),10,byid = TRUE)
##用它来剪切任何它横过
ii< - which(gIntersects(SLmap,SP,byid = TRUE))的SpatialLines段
#用分割版本替换违规行
# (但是当没有交叉点时跳过(例如,当lon_0 = 0时))
if(length(ii)){
SPii < - g Difference(SLmap [ii],SP,byid = TRUE)
SLmap< - rbind(SLmap [-ii],SPii)
}
return(SLmap)
}
## re-中心,并清理剩余的条纹
recenterAn dClean < - 函数(SLmap,lon_0){
llCRS <-CRS(+ proj = longlat + ellps = WGS84)## map包的CRS
newCRS < - CRS(paste( + proj = eqc + lon_0 =,lon_0,sep =))
##接收者
SLmap< - spTransform(SLmap,newCRS)
##识别剩余的'scratch-lines '通过搜索
##经过3条经线的2条线,间隔120度
v1 <-spTransform(readWKT(LINESTRING(-62 -89,-62 89)) (读取WKT(LINESTRING(58-89,58 89),p4s = 11CRS),newCRS)
v3< -spTransform(readWKT LINESTRING(178 -89,178 89),p4s = 11CRS),newCRS)
ii< - ((gIntersects(v1,SLmap,byid = TRUE)+
gIntersects(v2,SLmap ,byid = TRUE)+
gIntersects(v3,SLmap,byid = TRUE))> = 2)
SLmap [-ii]
}
##把它放在一起:
Recenter< - function(lon_0 = -100,grid = FALSE,...){
SLmap< - makeSLmap()
SLmap2< - sliceAtAntipodes(SLmap,lon_0)
recenterAndClean(SLmap2,lon_0)
}
##试试看
par(mfrow = c(2,2),mar = rep(1,4))
plot(Recenter(-90),col =grey40); box()##以90w
图为中心(Recenter(0),col =grey40); box()##以素数子午线
图为中心(Recenter(90),col =grey40); box()##以90e
图为中心(Recenter(180),col =grey40); box()##以国际日期行为中心
I am overlaying a world map from the
maps
package onto aggplot2
raster geometry. However, this raster is not centered on the prime meridian (0 deg), but on 180 deg (roughly the Bering Sea and the Pacific). The following code gets the map and recenters the map on 180 degree:require(maps) world_map = data.frame(map(plot=FALSE)[c("x","y")]) names(world_map) = c("lon","lat") world_map = within(world_map, { lon = ifelse(lon < 0, lon + 360, lon) }) ggplot(aes(x = lon, y = lat), data = world_map) + geom_path()
which yields the following output:
Quite obviously there are the lines draw between polygons that are on one end or the other of the prime meridian. My current solution is to replace points close to the prime meridian by NA, replacing the
within
call above by:world_map = within(world_map, { lon = ifelse(lon < 0, lon + 360, lon) lon = ifelse((lon < 1) | (lon > 359), NA, lon) }) ggplot(aes(x = lon, y = lat), data = world_map) + geom_path()
Which leads to the correct image. I now have a number of question:
- There must be a better way of centering the map on another meridian. I tried using the
orientation
parameter inmap
, but setting this toorientation = c(0,180,0)
did not yield the correct result, in fact it did not change anything to the result object (all.equal
yieldedTRUE
). - Getting rid of the horizontal stripes should be possible without deleting some of the polygons. It might be that solving point 1. also solves this point.
解决方案Here's a different approach. It works by:
- Converting the world map from the
maps
package into aSpatialLines
object with a geographical (lat-long) CRS. - Projecting the
SpatialLines
map into the Plate Carée (aka Equidistant Cylindrical) projection centered on the Prime Meridian. (This projection is very similar to a geographical mapping). - Cutting in two segments that would otherwise be clipped by left and right edges of the map. (This is done using topological functions from the
rgeos
package.) - Reprojecting to a Plate Carée projection centered on the desired meridian (
lon_0
in terminology taken from thePROJ_4
program used byspTransform()
in thergdal
package). - Identifying (and removing) any remaining 'streaks'. I automated this by searching for lines that cross g.e. two of three widely separated meridians. (This also uses topological functions from the
rgeos
package.)
This is obviously a lot of work, but leaves one with maps that are minimally truncated, and can be easily reprojected using
spTransform()
. To overlay these on top of raster images withbase
orlattice
graphics, I first reproject the rasters, also usingspTransform()
. If you need them, grid lines and labels can likewise be projected to match theSpatialLines
map.
library(sp) library(maps) library(maptools) ## map2SpatialLines(), pruneMap() library(rgdal) ## CRS(), spTransform() library(rgeos) ## readWKT(), gIntersects(), gBuffer(), gDifference() ## Convert a "maps" map to a "SpatialLines" map makeSLmap <- function() { llCRS <- CRS("+proj=longlat +ellps=WGS84") wrld <- map("world", interior = FALSE, plot=FALSE, xlim = c(-179, 179), ylim = c(-89, 89)) wrld_p <- pruneMap(wrld, xlim = c(-179, 179)) map2SpatialLines(wrld_p, proj4string = llCRS) } ## Clip SpatialLines neatly along the antipodal meridian sliceAtAntipodes <- function(SLmap, lon_0) { ## Preliminaries long_180 <- (lon_0 %% 360) - 180 llCRS <- CRS("+proj=longlat +ellps=WGS84") ## CRS of 'maps' objects eqcCRS <- CRS("+proj=eqc") ## Reproject the map into Equidistant Cylindrical/Plate Caree projection SLmap <- spTransform(SLmap, eqcCRS) ## Make a narrow SpatialPolygon along the meridian opposite lon_0 L <- Lines(Line(cbind(long_180, c(-89, 89))), ID="cutter") SL <- SpatialLines(list(L), proj4string = llCRS) SP <- gBuffer(spTransform(SL, eqcCRS), 10, byid = TRUE) ## Use it to clip any SpatialLines segments that it crosses ii <- which(gIntersects(SLmap, SP, byid=TRUE)) # Replace offending lines with split versions # (but skip when there are no intersections (as, e.g., when lon_0 = 0)) if(length(ii)) { SPii <- gDifference(SLmap[ii], SP, byid=TRUE) SLmap <- rbind(SLmap[-ii], SPii) } return(SLmap) } ## re-center, and clean up remaining streaks recenterAndClean <- function(SLmap, lon_0) { llCRS <- CRS("+proj=longlat +ellps=WGS84") ## map package's CRS newCRS <- CRS(paste("+proj=eqc +lon_0=", lon_0, sep="")) ## Recenter SLmap <- spTransform(SLmap, newCRS) ## identify remaining 'scratch-lines' by searching for lines that ## cross 2 of 3 lines of longitude, spaced 120 degrees apart v1 <-spTransform(readWKT("LINESTRING(-62 -89, -62 89)", p4s=llCRS), newCRS) v2 <-spTransform(readWKT("LINESTRING(58 -89, 58 89)", p4s=llCRS), newCRS) v3 <-spTransform(readWKT("LINESTRING(178 -89, 178 89)", p4s=llCRS), newCRS) ii <- which((gIntersects(v1, SLmap, byid=TRUE) + gIntersects(v2, SLmap, byid=TRUE) + gIntersects(v3, SLmap, byid=TRUE)) >= 2) SLmap[-ii] } ## Put it all together: Recenter <- function(lon_0 = -100, grid=FALSE, ...) { SLmap <- makeSLmap() SLmap2 <- sliceAtAntipodes(SLmap, lon_0) recenterAndClean(SLmap2, lon_0) } ## Try it out par(mfrow=c(2,2), mar=rep(1, 4)) plot(Recenter(-90), col="grey40"); box() ## Centered on 90w plot(Recenter(0), col="grey40"); box() ## Centered on prime meridian plot(Recenter(90), col="grey40"); box() ## Centered on 90e plot(Recenter(180), col="grey40"); box() ## Centered on International Date Line
这篇关于在绘制世界地图时使用与主要子午线不同的中心的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
- 将世界地图从