r - 如何在 map 上显示经度和纬度的值？

Question

我有一个投影为 EPSG:3410 的文件，我只想使用 R 绘制它。由于投影以米为单位，我对纬度和经度值(轴)有疑问.你可以在 map 上看到不是90 -90 180 -180而是其他数字。

绘制这张 map 的代码:

       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 替换 map 上的那些不寻常的值(因此原始文件没有变化)
• 如何去掉 map 上下的空间？

Answer 1

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

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))

要添加经纬度刻度，我们需要将 EPSG:3410(米)投影转换为经纬度(度)。我们可以在 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)

另一种选择是重新投影栅格，因此 map 是经纬度(度):

# 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)