在绘制世界地图时使用与主要子午线不同的中心

本文介绍了在绘制世界地图时使用与主要子午线不同的中心的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述 我将： lon = ifelse（lon lon = ifelse（（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投影，集中于所需子午线（程序中使用的术语中，code> lon_0 spTransform（）在 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 map2SpatialLines（wrld_p，proj4string = llCRS）} ## Clip SpatialLines整齐地沿着反对子午线 sliceAtAntipodes < - 函数（SLmap，lon_0）{ ##预备 long_180< - （lon_0 %% 360） - 180 llCRS< - CRS（+ proj = longlat + ellps = 地图对象的## CRS eqcCRS ##将地图重新​​映射到Equidistant圆柱/板块Caree投影 SLmap< ; - 线（cbind（long_180，c（-89,89））），ID（线）（cbind（long_180，c（-89,89）））; - 在子午线上形成一个窄的SpatialPolygon =cutter） SL SP ##用它来剪切任何它横过 ii< - which（gIntersects（SLmap，SP，byid = TRUE））的SpatialLines段＃用分割版本替换违规行＃ （但是当没有交叉点时跳过（例如，当lon_0 = 0时）） if（length（ii））{ SPii SLmap< - rbind（SLmap [-ii]，SPii）} return（SLmap）} ## re-中心，并清理剩余的条纹 recenterAn dClean llCRS newCRS ##接收者 SLmap< - spTransform（SLmap，newCRS） ##识别剩余的'scratch-lines '通过搜索 ##经过3条经线的2条线，间隔120度 v1 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 a ggplot2 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 in map, but setting this to orientation = c(0,180,0) did not yield the correct result, in fact it did not change anything to the result object (all.equal yielded TRUE).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 a SpatialLines 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 the PROJ_4 program used by spTransform() in the rgdal 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 with base or lattice graphics, I first reproject the rasters, also using spTransform(). If you need them, grid lines and labels can likewise be projected to match the SpatialLines 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" mapmakeSLmap <- 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 meridiansliceAtAntipodes <- 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 streaksrecenterAndClean <- 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 outpar(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 meridianplot(Recenter(90), col="grey40"); box() ## Centered on 90eplot(Recenter(180), col="grey40"); box() ## Centered on International Date Line 这篇关于在绘制世界地图时使用与主要子午线不同的中心的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！