问题描述
我想使用ggplot和LAmbert Azimuthal等面积投影（LAEA）映射一个水深的对象（来自“ marmap”包），以便我可以绘制按ID分组的点和路径。没关系，但是最终的地图投影变形了并且不适合LAEA。我尝试了不同的投影设置，但找不到问题（也许是在加固地图时？）。下面是我使用的代码。
首先，我加载地图并设置CRS：
library(marmap)
library(ggplot2)
library(dplyr)
# Get data -> raster world map (bathymetry map)
world <- getNOAA.bathy(-180,180,-90,90,resolution = 25,keep = TRUE)
# Switch to raster
world.ras <- as.raster(world)
# Set the projection to LAmbert Azimuthal Equal Area
my.proj <- "+proj=laea +lat_0=40 +lon_0=-30 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0 +units=m +no_defs"
world.ras.proj <- projectRaster(world.ras,crs = my.proj)
# Switch back to a bathy object
world.proj <- as.bathy(world.ras.proj)
# fortify raster map
world.proj.df <- fortify(world.proj)
带有位置数据的数据框如下所示：
hmm.stat <- data.frame(lat=c(77.16698,80.26016,59.98226,50.74124,62.19604,65.58466,80.03116,77.31931,51.29168,47.69987),lon=c(12.68404,27.43116,-37.875200,-48.6865620,-30.0061030,3.3246880,43.4332484,19.05524,-35.126510,-43.15675),ID=c(1,1,2,3,3))
coordinates(hmm.stat) <- c("lat","lon") #defined as coordinates
proj4string(hmm.stat) <- CRS("+init=epsg:4326") #defined CRS to WGS 84
hmm.stat.sp <- spTransform(hmm.stat,my.proj) #transform to LAEA
hmm.stat.df <- data.frame(hmm.stat.sp)
我在ggplot（）中使用的脚本给出了失真的投影：
map <- ggplot(world.proj.df,aes(x=x,y=y)) +
geom_raster(aes(fill=z),data=world.proj.df) +
geom_contour(aes(z=z),breaks=0,#contour for land
colour="black",size=1) +
scale_fill_gradientn(values = scales::rescale(c(-5000,3220)),colors = c("steelblue4","#C7E0FF","gray40","white")) +
coord_cartesian(ylim = c(1900000,9000000),xlim = c(150000,6600000),expand = FALSE) #set map extent
map + geom_path(data = hmm.stat.df,aes(x=lon,y=lat,group = ID),color ='gray10',alpha = 0.5,#transparency
size = 0.2) +
geom_point(data = hmm.stat.df,aes(x = lon,y = lat),colour = 'firebrick',size = 1.2,alpha = 0.3,#transparency
shape = 15) + #dots
ggtitle("Non-breeding movement") +
theme(panel.grid.minor.x=element_blank(),#remove grid
panel.grid.major.x=element_blank(),#remove grid
panel.grid.minor.y=element_blank(),#remove grid
panel.grid.major.y=element_blank(),#remove grid
axis.ticks.x = element_blank(),#remove axes graduation ticks
axis.ticks.y = element_blank(),#remove axes graduation ticks
panel.background = element_rect(fill = "white"),axis.text.x = element_blank(),axis.text.y = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank(),plot.title = element_text(colour = "grey25",size = 45,margin = margin(t = 0,r = 0,b = 12,l = 0)),panel.border = element_rect(colour = 'gray24',fill=NA,size=3),plot.margin = unit(c(0,0),"in"),# plot margins
legend.position="none") # remove all legends
我知道我可以直接使用plot（）而不是ggplot（），它似乎可以更好地处理栅格地图（请参见下文）。但是，我想使用ggplot（），因为它提供了更大的灵活性来按组处理地图和数据绘图。
使用plot（）进行映射而不会变形：
library(proj4)
library(graticule)
# Set colors for oceans and land masses
blues <- c("lightsteelblue4","lightsteelblue3","lightsteelblue2","lightsteelblue1")
greys <- c(grey(0.6),grey(0.93),grey(0.99))
# Projection of the coordinates of the corners of the area of interest
corners <- data.frame(lon = c(-50,-80,-10,-50),lat = c(5,5,5))
corners.proj <- project(as.matrix(corners),proj = my.proj)
par(mar = c(0,0))
plot(world.proj,image = TRUE,land = TRUE,lwd = 0.01,bpal = list(c(0,max(world.proj,na.rm = TRUE),greys),c(min(world.proj,blues)),deep = 0,shallow = 0,axes = FALSE,xlab = "",ylab = "",xlim = range(corners.proj[,1]),ylim = range(corners.proj[,2]))
# Add coastline
plot(world.proj,n = 1,lwd = 0.3,add = TRUE)
## specify where you want meridians and parallels
lons <- seq(-180,by = 20)
lats <- seq(0,80,by = 20)
## build the (projected graticules) lines
grat <- graticule(lons,lats,proj = my.proj)
## Add graticules to the plot
plot(grat,col = grey(0.93),add = TRUE,lty = 5)
非常感谢您的帮助！
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）