如何在地图上显示经纬度的值？

Question

我有一个投影为EPSG:3410 的文件，我只想使用 R 进行绘图。由于投影以米为单位，因此经纬度值（轴）存在问题。您可以在地图上看到除了 90 -90 180 -180 还有其他数字。

绘制此地图的代码：

       con1 <- file("C:\\Users\\data.bin","rb")
       r = raster(y)
      extent(r) = extent(c(xmn=-17334194,xmx=17334194,ymn=-7356860,ymx=7310585))
      plot(r)
      plot(wbuf, add=TRUE)

• 如何将地图上那些不寻常的值替换为 90 -90 180 -180（因此原始文件没有变化）
• 如何消除地图上下空间？

Answer 1

这里有一些解决方法可以避免空白并添加轴标签。由于两张地图都以米为单位，因此绘制经纬度标签并不容易。

library(raster) # for xmin, ymin etc.

# plot raster without axes and box
plot(r, asp=1, axes=F, box=F)

# add polygons
plot(wbuf, add=TRUE, axes=F)

# draw bounding box at raster limits 
rect(xmin(r), ymin(r), xmax(r), ymax(r))

要添加 lat-lon 刻度，我们需要将 EPSG:3410（米）投影转换为 lat-lon（度）。我们可以在 latlon 中创建虚拟空间点并将其转换为 3410，然后提取转换后的坐标以添加刻度

# y axis:
x = -180
y = seq(-90,90,30) # you can adapt this sequence if you want more or less labels
S <- SpatialPoints(cbind(x,y), proj4string = CRS("+proj=longlat +datum=WGS84"))
S2<- spTransform(S, CRS(" +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +no_defs ")) 
axis(side = 2, pos=xmin(r), at = S2@coords[,'y'], lab=y, las=1)


# x axis
y = -90
x = seq(-180,180,30)
S <- SpatialPoints(cbind(x,y), proj4string = CRS("+proj=longlat +datum=WGS84"))
S2<- spTransform(S, CRS(" +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +no_defs "))
axis(side = 1, pos=ymin(r), at = S2@coords[,'x'], lab=x)

另一种选择是重新投影栅格，因此地图是经纬度（度）：

# reproject and plot raster
r2 = projectRaster(r, crs=CRS("+proj=longlat +datum=WGS84"))
plot(r2, axes=F, box=F)

# there's an option to disable internal boundaries, which is nice
map(database = 'world', regions = '', interior = F, add=T)

# still have to add bounding box and axes
rect(xmin(r), ymin(r), xmax(r), ymax(r))
axis(1, pos = ymin(r2))
axis(2, pos= xmin(r2), las=1)