将空间坐标分配给R中的数组

Question

我从以下链接下载了一个数据文本文件：http://radon.unibas.ch/files/us_rn_50km.zip

解压缩后，我使用以下代码行绘制数据： # 加载库 库（字段）

# function to rotate a matrix (and transpose)
rotate <- function(x) t(apply(x, 2, rev))

# read data
data <- as.matrix(read.table("~/Downloads/us_rn_50km.txt", skip=6))
data[data<=0] <- NA

# rotate data
data <- rotate(data)

# plot data
mean.rn <- mean(data, na.rm=T)
image.plot(data, main=paste("Mean Rn emissions =", sprintf("%.3f", mean.rn)) )

这一切看起来都不错，但我希望能够在经纬度网格上绘制数据。我想我需要将此数组转换为 sp 类对象，但我不知道如何。我知道以下内容（来自网站）：“用于投影纬度和经度坐标的投影是用于北美地质十年 (DNAG) 地图的投影。投影类型是球形横轴墨卡托，基本纬度为零度和参考经度为 100 度 W。使用的比例因子为 0.926，没有假东移或北移。经纬度基准面是 NAD27，xy 坐标的单位是米。使用的椭球是 Clarke 1866。地图的分辨率为 50x50km”。但不知道如何处理这些数据。我试过了：

proj4string(data)=CRS("+init=epsg:4267")
data.sp <- SpatialPoints(data, CRS("+proj=longlat+ellps=clrk66+datum=NAD27") )

但有各种问题（与 NA 的问题），基本上我认为数据格式不正确。我认为 SpatialPoints 函数需要位置数据（二维）和第三个值数组与这些位置相关联（x、y、z 数据 - 我想我的问题是从我的数据中计算出 x 和 y！）

非常感谢任何帮助！

谢谢， 亚历克斯

Answer 1

有问题的文件是一个 ASCII 光栅网格。坐标在这种格式中是隐含的；标题描述（通常）左下角的位置，以及网格尺寸和分辨率。在此标题部分之后，由空格分隔的值描述了变量如何在网格中变化，值以行优先顺序给出。如果您有兴趣，请在文本编辑器中打开它。

您可以使用神奇的raster 包将此类文件导入 R，如下所示：

download.file('http://radon.unibas.ch/files/us_rn_50km.zip', 
          destfile={f <- tempfile()})
unzip(f, exdir=tempdir())
r <- raster(file.path(tempdir(), 'us_rn_50km.txt'))

您可以立即绘制它，而无需分配投影：

如果您不想将其转换为另一个 CRS，则不一定需要分配当前坐标系。但由于您确实想将其转换为 WGS84（地理），因此您需要先分配 CRS：

proj4string(r) <- '+proj=tmerc +lon_0=-100 +lat_0=0 +k_0=0.926 +datum=NAD27 +ellps=clrk66 +a=6378137 +b=6378137 +units=m +no_defs'

不幸的是，我不完全确定这个 proj4string 是否正确反映了website that provided the data 提供的信息（如果他们真的以标准格式提供了定义，那就太好了）。

分配CRS后，可以用projectRaster投影光栅：

r.wgs84 <- projectRaster(r, crs='+init=epsg:4326')

如果您愿意，可以将其写入您选择的光栅格式，例如：

writeRaster(r.wgs84, filename='whatever.tif')