如何查找经度/纬度网格线交点的坐标答案

【问题标题】：How to find coordinates of the longitude/latitude grid line intersection如何查找经度/纬度网格线交点的坐标
【发布时间】：2019-05-09 14:51:00
【问题描述】：

我想用本地地理投影绘制栅格，然后添加 lon/lat 网格。我需要在 xy 轴上的点上标记纬度/经度。让我展示一下方法：

library(raster)
linbrary(sp)
library(ggplot)
#local Albers geo projection 
aea <- CRS('+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 
             +x_0=4000000 +y_0=0 +datum=WGS84 +units=m +no_defs
             +ellps=WGS84 +towgs84=0,0,0 ')
wgs84 <- crs('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')

#create lon/lat grid spatiallinedataframe
lon_tick <- seq(-180,180,5)
lat_tick <- seq(-90,90,5)
j_line <- plyr::alply(lon_tick,1,function(x) cbind(x,lat_tick)) 
i_line <- plyr::alply(lat_tick,1,function(x) cbind(lon_tick,x) )
lines <- spLines(c(i_line,j_line),crs=wgs84)
lonlat_line <- lines %>% spTransform(.,aea_China)
lonlat_shp <- SpatialLinesDataFrame(lonlat_line, data = data.frame(ID=1))

#raster plot extent,xmin,xmax, ymin,ymax in order
ext <- c(3640676 ,5672676 ,3683208 ,5619208 )

#show the plot without raster tiles
p <- ggplot()+
 geom_polygon(data=lonlat_shp,aes(x=long,y=lat,group=group),
          fill=NA,color = "black",linetype=2) +
  coord_cartesian(xlim=ext[1:2],ylim=ext[3:4])

那么如何计算xy轴与lon/lat网格相交的点坐标（在albers投影下）？我认为在ggplot中使用scale_x_continuous时计算刻度坐标值后，标签xy轴很容易用lon/lat dgree打勾。

它喜欢求解方程：假设我想在 y 轴上获得一个坐标点 (x,Y)，与纬度 40°相交，x 在 ext 中称为 xmin，需要 Y。

aea_df <- data.frame( xx=x,yy=Y)
aea_point <- SpatialPoints(aea_df,aea)
lonlat_point <- spTransform(aea_point,wgs84) %>% as.data.frame

现在我们知道lonlat_point[1,2] 是 40 (40°N)，如何计算 Y？

【问题讨论】：

标签： r ggplot2

【解决方案1】：

我不知道 ggplot，但我认为这是获取坐标的方法。

示例数据

library(raster)
library(rgeos)
library(rgdal)
aea <- CRS('+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=4000000 +y_0=0 +datum=WGS84 +units=m +ellps=WGS84')
wgs84 <- crs('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
lonlat_lines <- as(as(raster(res=5), "SpatialPolygons"), "SpatialLines")
aea_lines <- spTransform(lonlat_lines, aea)
ext <- extent(c(3640676, 5672676, 3683208, 5619208))

解决方案

linext <- as(ext, "SpatialLines")
x <- gIntersection(aea_lines, linext)
xwgs <- spTransform(x, wgs84)
y <- coordinates(xwgs)

z <- crop(aea_lines, ext)
y <- coordinates(xwgs)
head(round(y,2))
#       x     y
#1 100.99 35.02
#1 100.72 40.01
#1 100.41 45.01
#1 100.06 50.00
#1 100.00 50.75
#1 100.00 51.77

现在进一步摆弄以获取它们各自轴的经度或纬度。

labs <- cbind(coordinates(x), y)
xlab <- labs[labs[,2] < 3740000, ]
ylab <- labs[labs[,1] < 3641000, ]
ylab <- ylab[abs(ylab[,4] - round(ylab[,4])) < .1, ]

还有情节

r <- as(raster(ext), "SpatialPolygons")
plot(r, col='light gray', border="white", xaxs="i", yaxs="i", las=1)
lines(z, lwd=2, lty=2)
axis(1, xlab[,1], xlab[,3], lwd=0, lwd.ticks=1)
axis(2, ylab[,2], round(ylab[,4]), lwd=0, lwd.ticks=1, las=1)

【讨论】：